Built with R 3.3.2

Bakgrund

En födelse-/dödsprocess1, (F/D), är en typ av Markovprocess definierad i kontinuerlig tid. Endast två typer av övergångar kan ske; “födelse”, vilket ökar tillståndsvariabeln med en enhet, resp “död” vilket på omvänt sätt minskar tillståndsvariabeln. Denna typ av process kan användas för att simulera olika naturfenomen, sjukdomar, bakterier och industriella processer och och andra kösystem. F/D-processens parametrar är “födelse-intensiteten” (ankomst) \(\displaystyle \lambda_{i}\), och “dödsintensiteten” (betjänings-intensiteten) \(\displaystyle {\mu _{i}}\). Enligt i uppgiften given förutsättning är båda dessa parametrar konstanta (tidsoberoende/stationära) för alla tidpunkter i.

I aktuell uppgift skall en F/D-process användas för att simulera en datorsystem, i det följande benämnt S. S består av 1 st CPU (“single server”) och 2 st seriekopplade buffertminnen. I dessa minnen kan inkommande arbeten kan lagras sekventiellt, i väntan på att få access till ledig CPU. Kön är således “trunkerad” till maximalt två arbeten. Tiden för inträffande av händelser: “ny uppgift”, är exponentialfördelad, \(\Delta_{idle} \sim Exp(1/\lambda)\). Betjäningstiden för servern är \(\Delta_{CPU}\sim Exp(1/\mu)\). S har således fyra tillstånd, utfallstrummet E = {0, 1, 2, 3}, där “0” betecknar att S är ledigt (idle), “1” att enbart CPU arbetar, “2” att CPU arbetar och en buffer är upptagen med ett väntande arbete, samt “3” att S är fullbelagt (båda buffrar är belagda med väntande arbeten och CPU arbetar). Aktuell typ av F/D-processen betecknas även M/M/1/2. Kön är i detta exempel begränsad till två platser, se Fig 0. Tillkommande, nya arbeten, inkommer enligt en Poissonprocess med intinsitet \(\lambda\). S har således en ändlig kö, med maximalt tre arbeten i väntan eller under processande. Ny arbeten ankomer oberoende av tidigare, kön har oberoende inkrement.

Fig 0. System S: M/M/1/2

Fig 0. System S: M/M/1/2

Sammanfattning

Teorin för Markovkedjor används för att lösa Laboration 3:s sex deluppgifter: 1. Systemet definieras 2. Kod för simulering av Makovkedja 3. Simulering av tre Förslag presenteras 4. Analys av simuleringsmodellens tider 5. Statistik över uppehållstid redovisas 6. Analytisk beräkning av ovanstående tider

Uppgift 1

Jämviktsfördelningen avseende sannnolikhet för tillståndsbyte definieras av en Poisson-process. M/M/1/2 är beteckningen för aktuellt system, bestående av en server (CPU) och ett tillhörande buffertminne, med buffertstorlek 2.

Systemets beskrivs av dess infinitesimala “generator”, G, även kallad Q-matris, där för varje rad balansekvationen ger att 0 = summa “slh-flöden”. Systemets intensitets-parametrar genererar element i en Q-matris, vilken för aktuell process blir:

\[ \mathbf{Q} : = \begin{pmatrix} -\lambda & \lambda& 0 & 0\\ \mu & -(\lambda + \mu) & \lambda & \ 0 \\ 0 & \mu & -(\lambda +\mu) & \lambda\\ 0 & 0 & \mu& -\mu \\ \end{pmatrix} \hskip 2cm \mathit{(i)} \] Från Q-matrisen erhålls övergångsmatrisen P, där varje rad har summan 1. P är en stokastisk matris:

\[ \mathbf{P} : = \begin{pmatrix} 0 & 1 & 0 & 0\\ \mu/(\lambda + \mu) & 0 & \lambda/(\lambda + \mu) & \ 0 \\ 0 & \mu/(\lambda + \mu) & 0 & \lambda/(\lambda + \mu)\\ 0 & 0 & 1 & 0 \\ \end{pmatrix} \hskip 2cm \mathit{(ii)} \]

Ankomsterna är enligt enligt givna förutsättningar Poissonfördelade. Tiden som åtgår för att i CPU processa inkommande arbeten är ochså Poissonfördela, med intenstet \(\mu\). Det innebär att tiden till nästa tillståndsbyte (upp eller ned), är exponentialfördelad (oberoende inkrement, minneslös), med intensitet \(\lambda +\mu\). Tiden till nästa tillståndsbyte, IA (InterArrivaltime), är då iid, d.v.s. identisk och likafördelad enligt \(\sim Exp(\lambda + \mu)\) \(\mathbf{ P(\mathrm{IA\leq t}) \ = 1- \exp^{-(\lambda + \mu)t}}.\) För grundtillståndet \(\mathrm E_{0}\) gäller dock att \(\mu = 0\), samt för \(\mathrm E_{3}\) att att \(\lambda = 0.\)

Sannolikheten att processen skall gå upp ett tillstånds-steg är lika med sannolikheten för en ny ankomst, minus sannolikheten för en sänkning av detsamma.

Sannolikheten för att lämna lägsta tillståndet (\(\mathrm E_{0}\), idle) under tiden t, : \(\mathbf P(E0 \rightarrow E1) = 1- \exp^{-\lambda t}\)

Sannolikheten att ytterligare öka tillståndet ett steg, från \(\mathrm E_{1}\) till \(\mathrm E_{2}\), respektive \(\mathrm E_{2}\) till \(\mathrm E_{3}\): \(\mathbf P\) (ökning, då i = 1 eller 2) = \(1- \exp^{-\lambda t}\)

Uppgift 2

Kod för att simulera en F/D-process, med egenskaper enligt uppgift 1

bd_process <- function(lambda, mu, initial_state = 0, steps = 500) {
  
# Simuleringen avbryts efter att ha nått antalet "steps" steg. *S* kan då befinna sig i vilket som helst av de 4 tillstånden. 
# Detta framgår t.ex. av Fig 3, där simuleringen avslutas då *S* befinner sig i "2", vilket medför att *S* stannar där och inte återgår till tillstånd "0".
# Intensitets-parametrarna: lambda_now, resp my_now, tilldelas aktuellt värde för vardera parametern. 
# Värdet kan vara lika med 0 vid MK start resp sluttillstånd. 
  time_now <- 0
  state_now <- initial_state

  # Dessa vektorer ska byggas på i loopen nedan
  time <- 0
  state <- initial_state

  for (i in 1:steps) {
    # Kontrollera systemets *S* aktuella tillstånd.
    # Om *S* är fullbelagt (state_now = 3) matas inte ny task in i *S* (lambda sätts till 0).
    # Annars ledigt är det tillgängligt för en ny task (med intensitet lambda)
    
    if (state_now == 3) {
        lambda_now <- 0
        } else {
          lambda_now <- lambda
        }
  # Om *S* är helt ledigt (idle) så finns kapacitet för nya tasks
  # my sätts  till 0 för (kan ej sänkas ytterligare från grundtillståndet)
  # Om inte det är fallet så avslutas pågående, (med intensitet my)
  # 
    if (state_now == 0) {
        mu_now <- 0
        } else {
          mu_now <- mu
    }
  #  Förfluten tid till nästa tillståndsändring 
  #  Tiden genereras genom att ge ett slumptal som argument till inversa fördelningen  F^-1 till ~(exp(my + lambda))
  time_to_transition <- -(1/(mu_now + lambda_now))*log(runif(1))

  # Undersök genererad sannolikhet för tillståndsändring: my * runif(1) > lambda * runif(1)
  # Om slh för avslut är större än slh för nystart går *S* in i nästa lägre tillstånd. 
  # I fall det omvända gäller går *S* sin in i nästa högre tillstånd.
  
  if (((mu_now + lambda_now )* runif(1)) < mu_now) {    # fler avslutade task än nytillkomna?
          state_now <- state_now - 1                    # Om ja: minska aktuellt tillstånd med ett
     } else {
          state_now <- state_now + 1                    # om nej: öka aktuellt tillstånd med ett
            }
      time_now <- time_now + time_to_transition # vektor med accumulerad tid mellan tillståndsändring,
                                                # med sista adderade element tidpunkten för nästa tillståndsändring.
      time <- c(time, time_now)                 # tidsvektorn, med sista element som är aktuell tid (tidpunkt för
                                                # sista tillståndsändringen)
      state <- c(state, state_now)              # Elementen i vektorn "state" består av benämningen (0, 1, 2, 3 =                                                                # tillståndsrummet) 
                                                # på de tillstånd som sekventiellt har genomlöpts, enl simulerad MK.
  }
  
# Returnera en lista med de två vektorerna tid och state
  
  list(tid = time, state = state, steps=steps)
}

Uppgift 3

Här nedan redovisas tre trappstegsdiagram, som återger tillståndsändringar vid simulering av F/D-processen enligt uppgift 2. Tre figurer presenteras, respektive figur representerar simulering med parametrar enligt Förslag 1, Förslag 2 och Förslag 3.

{set.seed(19500207)
  
#################################################
# Förslag 1 har processparameter, mu = 2, lambda = 10, enligt anvisning.
# Hämta ut vektorerna *tid* och *state*, samt antalet beräkningssteg steps
  
forslag1 <- bd_process(2, 10, 0, 100) # steps=100 ggr

time1 <- forslag1$tid
state1 <- forslag1$state
steps1 <- forslag1$steps


# Rita upp X(t) som funktion av t, 

plot(stepfun(time1[-1], state1),
  do.points = FALSE,
  xlab = "Tid (h)",
  ylab = "Tillstånd " ,
  main = expression(paste("Fig 1. Simulering av *S* enligt Förslag 1:  ", mu, "=10  ", lambda, "=2")),
  yaxt = "n",
  sub = paste("Tidssteg: från 0 till", steps1[1], " = antalet tillståndsändringar" ))
  axis(2, at = c(0, 1, 2, 3), las = 2)

######################################################

forslag2 <- bd_process(lambda = 6, mu = 10, 0, 100) # steps=100 ggr

time2 <- forslag2$tid
state2 <- forslag2$state
steps2 <- forslag2$steps

# Rita upp X(t) som funktion av t, 

plot(stepfun(time2[-1], state2),
  do.points = FALSE,
  xlab = "Tid (h)",
  ylab = "Tillstånd " ,
  main = expression(paste("Fig 2. Simulering av *S* enligt Förslag 2:  ", mu, "=10  ", lambda, "=6")),
  yaxt = "n",
  sub = paste("Tidssteg: från 0 till", steps2[1]))
  axis(2, at = c(0, 1, 2, 3), las = 2)

##################################################

forslag3 <- bd_process(lambda = 10, mu = 10, 0, 100) # steps=100 ggr

time3 <- forslag3$tid
state3 <- forslag3$state
steps3 <- forslag3$steps

# Rita upp X(t) som funktion av t, 

plot(stepfun(time3[-1], state3),
  do.points = FALSE,
  xlab = "Tid (h)",
  ylab = "Tillstånd " ,
  main = expression(paste("Fig 3. Simulering av *S* enligt Förslag 3:  ", mu, "=10  ", lambda, "=10")),
  yaxt = "n",
  sub = paste("Tidssteg: från 0 till", steps3[1]) )
  axis(2, at = c(0, 1, 2, 3), las = 2)
}

Uppgift 4

Nedan redovisas Tabell 1, som sammanställt återger “tiden” för att simulera 500 st tillståndsändringar för respektive förslag.

S_tabell tar sista elementet i vektorn time, med tidpunkter för tillståndsändringar. Detta sista element anger totaltiden för denna simulering.

Som framgår av tabell 1 så ökar S-tiden då aktiviteten minskar. Vid en relativt låg system-belastning sker endast få tillståndsändringar per tidsenhet. Det medför att S måste låtas arbeta motsvarande längre tid för att nå i Uppgift 4 ställt mål: 500 tillståndsändringar.

S_tabell <- data.frame(Antal_h_simulerad_tid = c(floor(time1[length(time1)]), trunc(time2[length(time2)]), round(time3[length(time3)], digits=0)))
rownames(S_tabell) <- c( "Förslag 1", "Förslag 2", "Förslag 3")
print(S_tabell)
          Antal_h_simulerad_tid
Förslag 1                    23
Förslag 2                     8
Förslag 3                     7
#S_tabell<- c(time1[length(time1)], time2[length(time2)], round(time3[length(time3)]))
Tabell 1: Antal S-timmar (h) för att simulera 500 tillståndsändringar (tidssteg)
Antal_h_simulerad_tid
Förslag 1 23
Förslag 2 8
Förslag 3 7

Uppgift 5

Andel av S uppehållstid i de fyra tillstånden analyseras. Förslag på lämpliga regler för nyttjande lämnas.

Del 1

Här redovisas koden, som för ett givet tillstånd, beräknar andelen av totala uppehållstiden (sojourntiden) som tillbringades i just detta tillstånd.

funktion proportion_in_state tar ett tillstånd s, valt ur utfallsrum (0, 1, 2 eller 3) och en simulerad enkel födelse/-dödsprocess, här specificerad genom function bd_process. function proportion_in_state returnerar andelen tid som processen spenderade i det angivna tillståndet s.

Del 2

För vart och ett av förslagen 1, 2, och 3, beräknas hur stor del av system-tiden som har tillbringats i respektive tillståndet. aatt är matrix (3x4), för att lagra andel av totala tiden i de 4 tillstånden för resp förslag.

            State0     State1      State2      State3
Förslag1 0.9082885 0.07928075 0.006131374 0.006299343
Förslag2 0.5018268 0.23489090 0.149956705 0.113325583
Förslag3 0.2496250 0.23610256 0.207728505 0.306543916

Tabell 2. Sammanställning av S kapacitetsutnyttjande för de tre förslagen. Simulering sker för 1000 tillståndsändringar (tidssteg).

df<- data.frame(Forslag = c(" 1: (lambda = 2, mu = 10)", " 2: (lambda = 6, mu = 10)", " 3: (lambda = 10, mu = 10)" ) )
df<- cbind(df, rbind(aatt[1,], aatt[2,], aatt[3,]))

names(df)[-1]<- paste0("Tillstånd ", 0:3)
knitr::kable(df, caption="Tabell 2: Andelen uppehållstid i de fyra (4) tilllstånden, redovisning sker för respektive förslag.", digits=4)
Tabell 2: Andelen uppehållstid i de fyra (4) tilllstånden, redovisning sker för respektive förslag.
Forslag Tillstånd 0 Tillstånd 1 Tillstånd 2 Tillstånd 3
1: (lambda = 2, mu = 10) 0.9083 0.0793 0.0061 0.0063
2: (lambda = 6, mu = 10) 0.5018 0.2349 0.1500 0.1133
3: (lambda = 10, mu = 10) 0.2496 0.2361 0.2077 0.3065

Kommentar till resultaten:

  • Förslag 1: innebär ett underutnyttjande av S-systemet. S har full beläggning endast ca 0,6 % av tiden det är i drift.

  • Förslag 2: Full beläggning sker ca 10 % av drift-tiden.

  • Förslag 3: innebär att S är fullt nyttjat ca 1/4 av den möjliga nuttjandetiden. S och ger då “felmeddelande” (var god vänta) för eventuella nya kunder. Det innebär en inte acceptabel servicenivå, långt lägre än det ställda kravt på högst 5 procent. (Detta är ett ofta förekommande resultat av “fritt tillträde” till en eftertraktad tjänst/vara/resurs, den s.k. “tradgey of commons”.)

Förslag 2 kan därför vara vara en lämplig “medelväg”, kombinerat med införande av prioritetsordning mellan de olika användarkategorierna, se nedan under rubrik Slutsats.

Uppgift 6

Här nedan anges en formel för den stationära fördelningen, som funktion av trafikintensiteten, \(\rho = \lambda / \mu\). Rekommendation lämnas beträffande val av vilken av de föreslagna reglerna som lämpligen bör införas.

Låt sannolikheten att befinna sig i tillstånd \(\mathrm E_{n}, \quad n =0, 1,... 4.\) betecknas \(\mathbf{P_{n}}\). Processen har parametern (intensitet) \(\lambda\) för att i följande tidssteg öka till nästa högre tillstånd och \(\mu\) är intensiteten för att minska till lägre tillstånd.

Under förutsättning att processen är “minneslös” fås följande balansekvationer:

\(\mathbf{P_{0}(t+ \Delta t) = (1-\lambda \Delta t)P_{0} + \mu \Delta P_{1} + o(\Delta t)}\)

\(\mathbf{P_{n}(t+ \Delta t) = \lambda \Delta t \, P_{n-1}(t) + (1 - (\lambda+ \mu))\Delta t \, P_{n}(t) + \mu \Delta t \,P_{n+1}(t) + o(\Delta t)}, \quad n= 1,2,3.\)

Vilket efter omskrivning ger följande uttryck: \[\mathbf { \frac{P_{n}(t+ \Delta t)- P_{n}(t)} {\Delta t} = \lambda P_{n-1}(t) -(\lambda + \mu)P_{n}(t) + \mu P_{n+1}(t) + o(\Delta t)/\Delta t} \quad (iii)\]

Efter gränsvärdesövergång, \(\lim{\Delta t \to 0}\), fås:

\(\mathbf{P_{n}^{´}(t) = \lambda P_{n-1}(t) - (\lambda + \mu) P_{n}(t) + \mu P_{n+1}(t)}, \quad n=1,2,3\)

Detta system av differential-differensekvationer benämns Kolmogorovs framåtekvationer2. Inom fysiken används den alternativa benämningen Fokker–Planck ekvationen. Ekvationssystemets begynnelsevärde ges av:

\(\mathbf{P_{0}^{´}(t) = - \lambda P_{0}(t) + \mu P_{1}(t)}\)

Enligt Kolmogorovs axiom3 gäller att summan av sannolikheter för tillståndsrummets samtliga tillstånd skall summera till 1.

Efter “lång” tid, när \(\lim{t \to \infty}\), har S intagit ett jämviktstillstånd, en jämviktsfördelning, då tidsderivatan (1/t),

\(\mathbf{P_{n}^{´}(t) = 0}\).

Därmed kan Kolmogorovs ekvationer förenklas enligt följande:

\(\mathbf {P_{1} = \frac{\lambda_{0}}{\mu_{1}}} \, \mathbf{ P_{0}}\)

\(\mathbf {P_{2} = \frac{\lambda_{1}}{\mu_{2}}} \, \mathbf {P_{1} = \frac{\lambda_{1}\lambda_{0}}{\mu_{2}\mu_{1}} } \, \mathbf{ P_{0}}\)

\(\mathbf P_{3} = \frac{\lambda_{2}}{\mu_{3}} \, \mathbf P_{2} = \frac{\lambda_{2}\lambda_{1}\lambda_{0}}{\mu_{3}\mu_{2}\mu_{1}} \, \mathbf P_{0} =\) \(\Pi_{i=1}^{3} (\lambda_{i-1}/\mu_{i}) \, \mathbf{ P_{0}} \hspace 2cm\)

Eftersom aktuell process är tidshomogen (stationär) kan förenkling göras, då fås: \(\lambda_{i} = \lambda, \; \mu_{i} =\mu \quad \forall\, i\), så erhålls

\(\mathbf {P_{0} = \frac {1} {1+\sum_{i=1}^{3} { (\lambda / \mu) } }} \hspace 7cm (iv)\)

Vi inför en beteckning för den s.k. trafikintensiteten, \(\rho = \lambda / \mu\), d.v.s kvoten av ankomst-/serviceintensiteten, fås uttrycket:

Det är inte nödvändigt att kräva (\(\rho\) < 1), eftersom kölängden är begränsad.
För att underlätta programmering (och slippa felmeddelande om resultatet NaN för rho =1) så behandlar vi fallet \(\rho = 1\) som ett särskilt fall. Då \(\rho = 1\), ser vi direkt att P1[4,] == 1/(N+1).

För varje Förslag (i) och tillstånd (j), applicerar vi formeln för en trunkerad geometric fördelning4 och sparar resutaten i matrisen P1:

\(\mathbf{ P1[i,j] = \frac{1-\rho} { 1-\rho^{N+1}} \rho^{j}}\)

Beräkning av värden till Tabell 3, enligt analytisk formel för jämviktsfördelningen:

N<- 3

P1<- matrix(nrow= 3, 
             ncol=4, 
             byrow <- TRUE)

# Förslag 1
rho <- 2/10
P1[1,1:4] <- rho^(0:3) * (1-rho)/(1-rho^(N+1))

# Förslag 2
rho <- 6/10
P1[2,1:4] <- rho^(0:3) * (1-rho)/(1-rho^(N+1))

Med nyttjande av ovan härled ekvation (iii), fås med aktuella parametrar följande sannolikheter:

För Förslag 1, med \(\lambda = 2\) och \(\mu = 10\) , fås trafikintensiteten \(\rho = \lambda / \mu =2/10 =1/5\):

\(\mathbf P_{0} = 1- \rho=\) 0.8012821,

\(\mathbf P_{1} = \rho P_{0} =\) 0.1602564,

\(\mathbf P_{2} = \rho^2 P_{0} =\) 0.0320513,

\(\mathbf P_{3} = \rho^3 P_{0} =\) 0.0064103.

Det vill säga en 0,64 procents risk för full beläggning av S.

För Förslag 2, med \(\lambda = 6\) och \(\mu = 10\) , fås \(\rho = \lambda / \mu =6/10 =3/5\):

\(\mathbf P_{0} = 1- \rho=\) 0.4595588,

\(\mathbf P_{1} = \rho P_{0} =\) 0.2757353,

\(\mathbf P_{2} = \rho^2 P_{0} =\) 0.1654412,

\(\mathbf P_{3} = \rho^3 P_{0} =\) 0.0992647.

Det vill säga 9,9 procents risk för full beläggning av S.

För Förslag 3, med \(\lambda = 10\) och \(\mu = 10\) , fås \(\rho = \lambda / \mu =10/10 = 1\):

\(\mathbf {P_{3}= P_{2}= P_{1}= P_{0} } = \rho^3/(1+3*\rho)= 1/4\).

Det betyder för S en 25 procents risk för överbeläggning (och därmed påföljande avvisning).

df<- data.frame(Forslag = c(" 1: (lambda = 2, mu = 10)", " 2: (lambda = 6, mu = 10)", " 3: (lambda = 10, mu = 10)" ) )
df<- cbind(df, rbind(P1[1,], P1[2,], P1[3,]))

names(df)[-1]<- paste0("Tillstånd ", 0:3)
knitr::kable(df, caption="Tabell 3: Stationär fördelning beräknad enligt formel, se Ross[^4]. Andelen uppehållstid i de fyra (4) tilllstånden, redovisas för respektive förslag.", digits=4)
Tabell 3: Stationär fördelning beräknad enligt formel, se Ross5. Andelen uppehållstid i de fyra (4) tilllstånden, redovisas för respektive förslag.
Forslag Tillstånd 0 Tillstånd 1 Tillstånd 2 Tillstånd 3
1: (lambda = 2, mu = 10) 0.8013 0.1603 0.0321 0.0064
2: (lambda = 6, mu = 10) 0.4596 0.2757 0.1654 0.0993
3: (lambda = 10, mu = 10) 1.0000 1.0000 1.0000 1.0000

Slutsats

Enligt Tabell 3 är risken för avvisande vid införande av Förslag 2 ca 10 procent. Detta är ett överskridande av det ställda kravet: högst 5 procent risk. Enligt gällande kravställning, maximalt 5 procent risk för avvisning, bör således Förslag 1 införas.

Enligt Tabell 3, Förslag 1, tar personalens beläggning dock endast en nästan försumbar del av systemets kapacitet: 0,64 procent av tiden är det risk att tillkommande arbeten avvisas.

Att avvisa samtliga studerande från den beräkningsresurs som S innebär, är dock knappast ett kostnadeffektivt nyttjande av systemet. Institutionsstyrelsen riskerar att få rykte som varande inkompetent och med bristande förmåga att administrera ett effektivt resursutnyttjande. Styrelsen bör därför överväga om tillträde och användning av S enligt ett reviderat Förslag 2, baserat på prioritetsordning mellan de två typerna av användare. Detta förslag, benämnt Förslag 4 innebär i så fall att endast personal tilldelas prioritet 1, och därmed alltid går före i arbetskön. Studenter bör ha tillgång till resursen, men med prioritet 2. Studenter bör kunna acceptera att personal går före i kön, utan att Studenterna känner sig kränkta. Prioritering av arbeten kan lösas tekniskt i buffert-systemet, genom att prioritet 1 arbeten alltid går in före, och vid behov “knuffar tillbaka” eller “knuffar ut” prio 2 arbeten från kön. I CPU pågående arbeten för prio 2 användare bör dock inte avbrytas, utan slutföras före prio 1 arbeten får överta CPU-kapaciteten. “Non-preemtive” prioritet förordas således före “absolut” prioritet (s.k. preemtive).

Referenser

  1. Ross, S.M. Introduction to Probability Models, s 359, 11 ed, 2014.↩

  2. Feller, W. An Introduction to Probability Theory and Its Applications, s 407, Vol 1, 2 ed, 1957.↩

  3. Kolmogorov A. N. Foundations of the Theory of Probability, s 2, 2 ed, 1956.↩

  4. Privault N. Understanding Markov Chains, s 195, 2013.↩

  5. Ross, S.M. ibid, s 496.↩