Zacznijmy od najprostszej sytuacji. Rozważmy skończony zbiór X, który będziemy nazywali przestrzenią stanów. Czasem wygodnie przyjąć, że X=1,2,…,d. Jest to tylko umowne ponumerowanie stanów.

Jeśli Xn=x, to mówimy, że łańcuch w chwili n znajduje się w stanie x∈X. Warunek (i) w Definicji 7.1 znaczy tyle, że przyszła ewolucja łańcucha zależy od stanu obecnego, ale nie zależy od przeszłości. Łańcuch jest jednorodny, jeśli prawo ewolucji łańcucha nie zmienia się w czasie. W dalszym ciągu rozpatrywać będziemy głównie łańcuchy jednorodne. Macierz

P=P⁢x,yx,y∈X=P⁢1,1⋯P⁢1,y⋯P⁢1,d⋮⋱⋮⋱⋮P⁢x,1⋯P⁢x,y⋯P⁢x,d⋮⋱⋮⋱⋮P⁢d,1⋯P⁢d,y⋯P⁢d,d

nazywamy macierzą (prawdopodobieństw) przejścia łańcucha. Jest to macierz stochastyczna, to znaczy P⁢x,y≥0 dla dowolnych stanów x,y∈X oraz ∑yP⁢x,y=1 dla każdego x∈X. Jeśli q⁢x=PX0=x, to wektor wierszowy

q⊤=q⁢1,…,q⁢x,…,q⁢d

nazywamy rozkładem początkowym łańcucha (oczywiście, ∑xq⁢x=1). Jest jasne, że

P(X0=x0,X1=x1,…,Xn-1=xn-1,Xn=xn)=q(x0)P(x0,x1)⋯P(xn-1,xn).

(7.1)

Ten wzór określa jednoznacznie łączny rozkład prawdopodobieństwa zmiennych X0,X1,…,Xn,… i może być przyjęty za definicję jednorodnego łańcucha Markowa. W tym sensie możemy utożsamić łańcuch z parą q,P: łączny rozkład prawdopodobieństwa jest wyznaczony przez podanie rozkładu początkowego i macierzy przejścia.

Opiszemy teraz bardzo ogólną konstrukcję, która jest podstawą algorytmów generujących łańcuchy Markowa. Wyobraźmy sobie, jak zwykle, że mamy do dyspozycji ciąg U=U0,U1,…,Un,… ,,liczb losowych”, produkowanych przez komputerowy generator, czyli z teoretycznego punktu widzenia ciąg niezależnych zmiennych losowych o jednakowym rozkładzie, U⁢0,1. Niech ϕ:0,1→X i ψ:X×0,1→X będą takimi funkcjami, że

⁢⁢Pψ⁢U=x=q⁢x⁢⁢dla każdego ⁢x∈X,⁢⁢Pϕ⁢x,U=y=P⁢x,y⁢⁢dla dowolnych ⁢x,y∈X.⁢

(7.2)

Jeśli określimy X0,X1,…,Xn,… rekurencyjnie w następujący sposób:

X0=ψ⁢U0,⁢Xn+1=ϕ⁢Xn,Un+1,

(7.3)

to jest jasne, że tak otrzymany ciąg zmiennych losowych jest łańcuchem Markowa z rozkładem początkowym q i macierzą przejścia P.

Na zakończenie tego podrozdziału zróbmy kilka prostych uwag.

• Wszystko, co w tym podrozdziale zostało powiedziane można niemal mechanicznie uogólnić na przypadek nieskończonej ale przeliczalnej przestrzeni X.

• Można sobie wyobrażać, że dla ustalonego x, zbiór u:ϕ⁢x,u=y jest odcinkiem długości P⁢x,y, ale nie musimy tego żądać. Nie jest istotne, że zmienne Ui mają rozkład U⁢0,1. Moglibyśmy założyć, że są określone na dość dowolnej przestrzeni U. Istotne jest, że te zmienne są niezależne, mają jednakowy rozkład i zachodzi wzór (7.2). W praktyce bardzo często do zrealizowania jednego ,,kroku” łańcucha Markowa generujemy więcej niż jedną ,,liczbę losową”.

Przypadek ciągłej lub raczej ogólnej przestrzeni stanów jest w istocie równie prosty, tylko oznaczenia trzeba trochę zmienić i sumy zamienić na całki. Dla prostoty przyjmijmy, że X⊆Rd. Poniżej sformułujemy definicję w taki sposób, żeby podkreślić analogię do przypadku dyskretnego i pominiemy subtelności teoretyczne. Zwróćmy tylko uwagę, że prawdopodobieństwo warunkowe nie może tu być zdefiniowane tak elementarnie jak w Definicji 7.1, bo zdarzenie warunkujące może mieć prawdopodobieństwo zero.

Funkcja P(⋅,⋅), której argumentami są punkt x i zbiór B, jest nazywana jądrem przejścia. Ważne dla nas jest to, że dla ustalonego x∈X, jądro P(x,⋅) rozważane jako funkcja zbioru jest rozkładem prawdopodobieństwa. W wielu przykładach jest to rozkład zadany przez gęstość,

P⁢x,B=∫Bp⁢x,y⁢d⁢y.

Wtedy odpowiednikiem wzoru (7.1) jest następujący wzór na łączną gęstość:

⁢p⁢x0,x1,…,xn-1,xn=q⁢x0⁢p⁢x0,x1⁢⋯⁢p⁢xn-1,xn,

(7.4)

gdzie q(⋅) jest gęstością rozkładu początkowego, zaś p(xi,xi+1)=p(xi+1|xi) jest gęstością warunkową. Sformułowaliśmy Definicję 7.2 nieco ogólniej głównie dlatego, że w symulacjach ważną rolę odgrywają łańcuchy dla których prawdopodobieństwo przejścia nie ma gęstości. Przykładem może być łańcuch, który z niezerowym prawdopodobieństwem potrafi ,,stać w miejscu”, to znaczy P(x,{x})=:α(x)>0, i ma prawdopodobieństwo przejścia postaci, powiedzmy

P(x,B)=α(x)I(x∈B)+(1-α(x))∫Bp(x,y)dy.

Łańcuch Markowa na ogólnej przestrzeni stanów można rekurencyjnie generować w sposób opisany wzorem (7.3), to znaczy X0=ψ⁢U0 i Xn+1=ϕ⁢Xn,Un+1 dla ciągu i.i.d. U0,U1,…,Un,… z rozkładu U⁢0,1. Niech Q(⋅) oznacza rozkład początkowy. Funkcje ψ i ϕ muszą spełniać warunki

⁢⁢⁢Pψ⁢U∈B=Q⁢B⁢⁢dla każdego (mierzalnego) ⁢B⊆X,⁢⁢Pϕ⁢x,U∈B=P⁢x,B⁢⁢dla dowolnych ⁢x∈X,B⊆X.⁢

(7.5)

Różnica między (7.2) i (7.5) jest tylko taka, że w ogólnej sytuacji rozkłady prawdopodobieństwa zmiennych losowych ψ⁢U i ϕ⁢x,U, czyli odpowiednio Q(⋅) i P(x,⋅) nie muszą być dyskretne. Ale można zastosować w zasadzie każdą metodę generacji zmiennych o zadanym rozkładzie, przy czym ten zadany rozkład zależy od x.

Rozważmy proces stochastyczny X⁢t,t⩾0, czyli rodzinę zmiennych losowych o wartościach w skończonej przestrzeni stanów X, indeksowaną przez parameter t, który interpretujemy jako ciągły czas. Pojedyncze stany oznaczamy przez x,y,z,…∈X lub podobnie. Załóżmy, że trajektorie są prawostronnie ciągłymi funkcjami mającymi lewostronne granice (prawie na pewno).

Sens założenia (i) jest taki sam jak w Definicji 7.1. Żądamy mianowicie, żeby zmienna X⁢t+s byłą warunowo niezależna od zmiennych X⁢u,0⩽u<s pod warunkiem zdarzenia X⁢s=x. Oczywiście, Pt⁢x,y jest prawdopodobieństwem przejścia w ciągu czasu t. Rozkład początkowy oznaczymy przez q⁢x=PX⁢0=x. Ponieważ X jest skończona, q może być przedstawiony jako wektor a Pt jako macierz.

Odpowiednikiem macierzy przejścia jest macierz intensywności przejść zdefiniowana następująco.

Q⁢x,y=limt→0⁡1h⁢Ph⁢x,y-I⁢x,y,

gdzie I=P0 jest macierzą identycznościową. Można udowodnić, że granice istnieją. Tak więc Q=dd⁢t⁢Ptt=0 i Ph=I+h⁢Q+o⁢h przy h↘0. Innymi słowy, Q⁢x,y jest rzeczywiście ,,intensywnością skoków z x do y” (na jednostkę czasu). Mamy wzór analogiczny do (6.3):

P(X(t+h)=y|X(t)=x)=hQ(x,y)+o(h) dla x≠y;P(X(t+h)=x|X(t)=x)=1+hQ(x,x)+o(h),h↘0.

Zauważmy, że mamy ∑yQ⁢x,y=0. Wspomnijmy jeszcze mimochodem, że macierze przejścia wyrażają się przez macierz intensywności Q przy pomocy ,,macierzowej funkcji wykładniczej”: Pt=exp⁡t⁢Q. Dla uproszczenia notacji, niech Q⁢x=-Q⁢x,x=∑y≠xQ⁢x,y oznacza “intensywność wszystkich skoków ze stanu x”.

Z tego co zostało powiedziane łatwo się domyślić, że generowanie procesu Markowa można zorganizować podobnie jak dla procesu Poissona, poprzez czasy skoków. Niech 0<T1<T2<⋯<Tn<⋯ będą kolejnymi momentami skoków,

⁢T1⁢=inf⁡t>0:X⁢t≠X⁢0;⁢Tn+1⁢=inf⁡t>Tn:X⁢t≠X⁢Tn.

Jeśli X⁢Tn=x, to czas oczekiwania na następny skok, Wn+1=Tn+1-Tn zależy tylko od x i ma rozkład wykładniczy z parametrem Q⁢x. W szczególności, mamy E(Wn+1|X(Tn)=x)=1/Q(x). Jeśli obserwujemy proces w momentach skoków, czyli skupimy uwagę na ciągu X⌃n=X⁢Tn to otrzymujemy łańcuch Markowa z czasem dyskretnym, nazywany czasem ,,szkieletem”. Jego prawdopodobieństwa przejścia są takie:

P⌃(x,y)=P(X⌃n+1=y|X⌃n=x)=Q⁢x,y/Q⁢xjeśli x≠y;0jeśli x=y.

(7.6)

Oczywiście, cała trajektoria procesu X⁢t jest w pełni wyznaczona przez momenty skoków T1,…,Tn,… i kolejne stany łańcucha szkieletowego chain X⁢0=X⌃0,X⌃1,…,X⌃n,…. Po prostu, X⁢t=X⌃n dla Tn≤t<Tn+1. Następujący algorytm formalizuje opisany powyżej sposób generacji.

Rzecz jasna, w praktyce trzeba ustalić sobie skończony horyzont czasowy h i przerwać symulacje gdy Ti>h. Zauważmy, że ten algorytm bez żadnych modyfikacji nadaje się do symulowania procesu Markowa na nieskończoniej ale przeliczalnej przestrzeni stanów, na przykład X=0,1,2,….

Czasami warto zmodyfikować algorytm w następujący sposób. Generowanie procesu Markowa można rozbić na dwie fazy: najpierw wygenerować ,,potencjalne czasy skoków” a następnie symulować ,,szkielet”, który będzie skakał wyłącznie w poprzednio otrzymanych momentach (ale nie koniecznie we wszystkich spośród nich). Opiszemy dokładniej tę konstrukcję, która bazuje na własnościach procesu Poissona. Wyobraźmy sobie najpierw, że dla każdego stanu x∈X symulujemy niezależnie jednorodny proces Poissona Rx o punktach skoku R1x<⋯<Rkx<⋯ przy czym ten proces ma intensywność Q⁢x. Jeśli teraz, w drugiej fazie X⌃i-1=X⁢Ti-1=x, to za następny moment skoku wybierzemy najbliższy punkt procesu Rx na prawo od Ti-1, czyli Ti=min⁡Rkx:Rkx>Ti-1. Z własności procesu Poissona (patrz Stwierdzenie 6.2 i jego dowód) wynika, że zmienna Ti-Ti-1 ma rozkład wykładniczy z parametrem Q⁢x i metoda jest poprawna. Naprawdę nie ma nawet potrzeby generowania wszystkich procesów Poissona Rx. Warto posłużyć się metodą przerzedzania i najpierw wygenerować jeden proces o dostatecznie dużej intensywności a następnie ,,w miarę potrzeby” go przerzedzać. Niech R* będzie procesem Poissona o intensywności Q*⩾maxx⁡Q⁢x i oznaczmy jego punkty skoków przez R1<⋯<Rk<⋯. Jeśli w drugiej fazie algorytmu mamy X⁢Ri-1=x w pewnym momencie Ri-1∈R*, to zaczynamy przerzedzać proces R* w taki sposób, aby otrzymać proces o intensywności Q⁢x⩽Q*. Dla każdego punktu Rj, j⩾i powinniśmy rzucić monetą i z prawdopodobieństwem 1-Q⁢x/Q* ten punkt usunąć. Ale jeśli usuniemy punkt Ri to znaczy, że w tym punkcie nie ma skoku, czyli X⁢Ri=x. Jeśli punkt Ri zostawimy, to wykonujemy skok, losując kolejny stan z prawdopodobieństwem P⌃(x,⋅). W rezultacie, stan X⁢Ri jest wylosowany z rozkładu prawdopodobieństwa

⁢P~⁢x,y=Q⁢x,y/Q*jeśli x≠y;1-Q⁢x/Q*jeśli x=y.

(7.7)

Niech X~n=X⁢Ri będzie ,,nadmiarowym szkieletem” procesu. Jest to łańcuch Markowa, który (w odróżnieniu od ,,cieńszego szkieletu” X⌃n) może w jednym kroku pozostać w tym samym stanie i ma prawdopodobieństwa przejścia P(X~n=y|X~n-1=x)=P~(x,y).

Odpowiada temu następujący dwufazowy algorytm.

• Faza I

• Faza II

Algorytm jest bardzo podobny do metody wynalezionej przez Gillespie'go do symulowania wielowymiarowych procesów urodzin i śmierci (głównie w zastosowaniach do chemii). Te procesy mają najczęściej przestrzeń stanów postaci X=0,1,2,…d; stanem układu jest układ liczb x=x⁢1,…,x⁢d, gdzie x⁢i jest liczbą osobników (na przykład cząstek) typu i. Przestrzeń jest nieskończona, ale większa część rozważań przenosi się bez trudu. Pojawia się tylko problem w ostatnim algorytmie jeśli maxx⁡Q⁢x=∞, czego nie można wykluczyć. Zeby sobie z tym poradzić, można wziąć Q* dostatecznie duże, żeby ,,na ogół” wystarczało, a w mało prawdopodobnym przypadku zawitania do stanu x z Q⁢x>Q* przerzucać się na pierwszy, jednofazowy algorym.