Zagadnienia

7.1 Czas dyskretny, przestrzeń dyskretna
7.2 Czas dyskretny, przestrzeń ciągła
7.3 Czas ciągly, przestrzeń dyskretna

7. Symulowanie procesów stochastycznych II. Procesy Markowa

7.1. Czas dyskretny, przestrzeń dyskretna

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

Definicja 7.1

(i) Ciąg $X_{0},X_{1},\ldots,X_{n},\ldots$ zmiennych losowych o wartościach w ${\cal X}$ nazywamy łańcuchem Markowa, jeśli dla każdego $n=1,2,\ldots$ i dla każdego ciągu $x_{0},x_{1},\ldots,x_{n},x_{{n+1}}$ punktów przestrzeni ${\cal X}$ ,

$\begin{split}&\mathbb{P}(X_{{n+1}}=x_{{n+1}}|X_{n}=x_{n},X_{{n-1}}=x_{{n-1}},\ldots,X_{1}=x_{1},X_{0}=x_{0})\\ &=\mathbb{P}(X_{{n+1}}=x_{{n+1}}|X_{n}=x_{n}),\\ \end{split}$

(o ile tylko $\mathbb{P}(X_{n}=x_{n},X_{{n-1}}=x_{{n-1}},\ldots,X_{1}=x_{1},X_{0}=x_{0})>0$ , czyli prawdopodobieństwo warunkowe w tym wzorze ma sens).

(ii) Łańcuch Markowa nazywamy jednorodnym, jeśli dla dowolnych stanów $x$ i $y$ i każdego $n$ możemy napisać

$\mathbb{P}(X_{{n+1}}=y|X_{n}=x)=P(x,y),$

to znaczy prawdopodobieństwo warunkowe w powyższym wzorze zależy tylko od $x$ i $y$ , ale nie zależy od $n$ .

Jeśli $X_{n}=x$ , to mówimy, że łańcuch w chwili $n$ znajduje się w stanie $x\in{\cal 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=\left(P(x,y)\right)_{{x,y\in{\cal X}}}=\left(\begin{array}[]{ccccc}P(1,1)&\cdots&P(1,y)&\cdots&P(1,d)\cr\vdots&\ddots&\vdots&\ddots&\vdots\cr P(x,1)&\cdots&P(x,y)&\cdots&P(x,d)\cr\vdots&\ddots&\vdots&\ddots&\vdots\cr P(d,1)&\cdots&P(d,y)&\cdots&P(d,d)\cr\end{array}\right)$

nazywamy macierzą (prawdopodobieństw) przejścia łańcucha. Jest to macierz stochastyczna, to znaczy $P(x,y)\geq 0$ dla dowolnych stanów $x,y\in{\cal X}$ oraz $\sum _{y}P(x,y)=1$ dla każdego $x\in{\cal X}$ . Jeśli $q(x)=\mathbb{P}(X_{0}=x)$ , to wektor wierszowy

$q^{\top}=\big(q(1),\ldots,q(x),\ldots,q(d)\big)$

nazywamy rozkładem początkowym łańcucha (oczywiście, $\sum _{x}q(x)=1$ ). Jest jasne, że

$\begin{split}&\mathbb{P}(X_{0}=x_{0},X_{1}=x_{1},\ldots,X_{{n-1}}=x_{{n-1}},X_{n}=x_{n})\\ &=q(x_{0})P(x_{0},x_{1})\cdots P(x_{{n-1}},x_{n}).\end{split}$

(7.1)

Ten wzór określa jednoznacznie łączny rozkład prawdopodobieństwa zmiennych $X_{0},X_{1},\ldots,X_{n},\ldots$ 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=U_{0},U_{1},\ldots,U_{n},\ldots$ ,,liczb losowych”, produkowanych przez komputerowy generator, czyli z teoretycznego punktu widzenia ciąg niezależnych zmiennych losowych o jednakowym rozkładzie, ${\rm U}(0,1)$ . Niech $\phi:[0,1]\to{\cal X}$ i $\psi:{\cal X}\times[0,1]\to{\cal X}$ będą takimi funkcjami, że

$\begin{split}&\mathbb{P}\left(\psi(U)=x\right)=q(x)\quad\textrm{dla każdego }x\in{\cal X},\\ &\mathbb{P}\left(\phi(x,U)=y\right)=P(x,y)\quad\textrm{dla dowolnych }x,y\in{\cal X}.\\ \end{split}$

(7.2)

Jeśli określimy $X_{0},X_{1},\ldots,X_{n},\ldots$ rekurencyjnie w następujący sposób:

$X_{0}=\psi(U_{0}),\qquad X_{{n+1}}=\phi\left(X_{n},U_{{n+1}}\right),$

(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 ${\cal X}$ .
Można sobie wyobrażać, że dla ustalonego $x$ , zbiór $\{ u:\phi(x,u)=y\}$ jest odcinkiem długości $P(x,y)$ , ale nie musimy tego żądać. Nie jest istotne, że zmienne $U_{{i}}$ mają rozkład ${\rm U}(0,1)$ . Moglibyśmy założyć, że są określone na dość dowolnej przestrzeni ${\cal 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ą”.

7.2. Czas dyskretny, przestrzeń ciągła

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 ${\cal X}\subseteq\mathbb{R}^{{d}}$ . 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.

Definicja 7.2

(i) Ciąg $X_{0},X_{1},\ldots,X_{n},\ldots$ zmiennych losowych o wartościach w ${\cal X}$ nazywamy łańcuchem Markowa, jeśli dla każdego $n=1,2,\ldots$ , dla każdego ciągu $x_{0},x_{1},\ldots,x_{n}$ punktów przestrzeni ${\cal X}$ oraz dowolnego (borelowskiego ) zbioru $B\subseteq{\cal X}$ ,

$\begin{split}&\mathbb{P}(X_{{n+1}}\in B|X_{n}=x_{n},X_{{n-1}}=x_{{n-1}},\ldots,X_{1}=x_{1},X_{0}=x_{0})\\ &=\mathbb{P}(X_{{n+1}}\in B|X_{n}=x_{n}),\\ \end{split}$

(dla prawie wszystkich punktów $(x_{n},x_{{n-1}},\ldots,x_{1},x_{0})$ ).

(ii) Łańcuch Markowa nazywamy jednorodnym, jeśli dla dowolnego zbioru $B\subseteq{\cal X}$ i każdego $n$ możemy napisać

$\mathbb{P}(X_{{n+1}}\in B|X_{n}=x)=P(x,B)$

(dla prawie wszystkich punktów $x$ ).

Funkcja $P(\cdot,\cdot)$ , 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\in{\cal X}$ , jądro $P(x,\cdot)$ 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)=\int _{B}p(x,y){\rm d}y.$

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

$\qquad p(x_{0},x_{1},\ldots,x_{{n-1}},x_{n})=q(x_{0})p(x_{0},x_{1})\cdots p(x_{{n-1}},x_{n}),$

(7.4)

gdzie $q(\cdot)$ jest gęstością rozkładu początkowego, zaś $p(x_{{i}},x_{{i+1}})=p(x_{{i+1}}|x_{{i}})$ 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\})=:\alpha(x)>0$ , i ma prawdopodobieństwo przejścia postaci, powiedzmy

$P(x,B)=\alpha(x)\mathbb{I}(x\in B)+(1-\alpha(x))\int _{B}p(x,y){\rm d}y.$

Łańcuch Markowa na ogólnej przestrzeni stanów można rekurencyjnie generować w sposób opisany wzorem (7.3), to znaczy $X_{0}=\psi(U_{0})$ i $X_{{n+1}}=\phi\left(X_{n},U_{{n+1}}\right)$ dla ciągu i.i.d. $U_{{0}},U_{{1}},\ldots,U_{{n}},\ldots$ z rozkładu ${\rm U}(0,1)$ . Niech $Q(\cdot)$ oznacza rozkład początkowy. Funkcje $\psi$ i $\phi$ muszą spełniać warunki

$\qquad\begin{split}&\mathbb{P}\left(\psi(U)\in B\right)=Q(B)\quad\textrm{dla każdego (mierzalnego) }B\subseteq{\cal X},\\ &\mathbb{P}\left(\phi(x,U)\in B\right)=P(x,B)\quad\textrm{dla dowolnych }x\in{\cal X},B\subseteq{\cal X}.\\ \end{split}$

(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 $\psi(U)$ i $\phi(x,U)$ , czyli odpowiednio $Q(\cdot)$ i $P(x,\cdot)$ 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$ .

7.3. Czas ciągly, przestrzeń dyskretna

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

Definicja 7.3

(i) Mówimy, że $X(t)$ jest procesem Markowa, jeśli dla dowolnych $x,y\in{\cal X}$ oraz $s,t\geqslant 0$ ,

$\mathbb{P}\left(X(s+t)=y|X(s)=x,X(u),0\leqslant u<s\right)=\mathbb{P}\left(X(t+s)=y|X(s)=x\right).$

(ii) Proces Markowa nazywamy jednorodnym, jeśli dla dowolnych stanów $x$ i $y$ i każdego $t$ i $s$ możemy napisać

$\mathbb{P}\left(X(s+t)=y|X(s)=x\right)=P^{t}(x,y),$

to znaczy prawdopodobieństwo warunkowe w powyższym wzorze zależy tylko od $x$ , $y$ i $t$ , ale nie zależy od $s$ .

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\leqslant u<s$ pod warunkiem zdarzenia $X(s)=x$ . Oczywiście, $P^{t}(x,y)$ jest prawdopodobieństwem przejścia w ciągu czasu $t$ . Rozkład początkowy oznaczymy przez $q(x)=\mathbb{P}\left(X(0)=x\right)$ . Ponieważ ${\cal X}$ jest skończona, $q$ może być przedstawiony jako wektor a $P^{t}$ jako macierz.

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

$Q(x,y)=\lim _{{t\to 0}}\frac{1}{h}\left[P^{h}(x,y)-I(x,y)\right],$

gdzie $I=P^{0}$ jest macierzą identycznościową. Można udowodnić, że granice istnieją. Tak więc $Q=\frac{{\rm d}}{{\rm d}t}P^{t}|_{{t=0}}$ i $P^{h}=I+hQ+o(h)$ przy $h\searrow 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):

$\begin{split}\mathbb{P}\left(X(t+h)=y|X(t)=x\right)&=hQ(x,y)+o(h)\text{ dla $x\not=y$};\\ \mathbb{P}\left(X(t+h)=x|X(t)=x\right)&=1+hQ(x,x)+o(h),\quad h\searrow 0.\end{split}$

Zauważmy, że mamy $\sum _{y}Q(x,y)=0$ . Wspomnijmy jeszcze mimochodem, że macierze przejścia wyrażają się przez macierz intensywności $Q$ przy pomocy ,,macierzowej funkcji wykładniczej”: $P^{t}=\exp[tQ]$ . Dla uproszczenia notacji, niech $Q(x)=-Q(x,x)=\sum _{{y\not=x}}Q(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<T_{1}<T_{2}<\cdots<T_{n}<\cdots$ będą kolejnymi momentami skoków,

$\begin{split} T_{1}&=\inf\left\{ t>0:X(t)\not=X(0)\right\};\\ T_{{n+1}}&=\inf\left\{ t>T_{n}:X(t)\not=X(T_{n})\right\}.\end{split}$

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

$\qquad{\hat{P}}(x,y)=\mathbb{P}({\hat{X}}_{{n+1}}=y|{\hat{X}}_{{n}}=x)=\begin{cases}Q(x,y)/Q(x)&\text{jeśli $x\not=y$;}\\ 0&\text{jeśli $x=y$}.\end{cases}$

(7.6)

Oczywiście, cała trajektoria procesu $X(t)$ jest w pełni wyznaczona przez momenty skoków $T_{1},\ldots,T_{n},\dots$ i kolejne stany łańcucha szkieletowego chain $X(0)={\hat{X}}_{{0}},{\hat{X}}_{{1}},\ldots,{\hat{X}}_{{n}},\ldots$ . Po prostu, $X(t)={\hat{X}}_{{n}}$ dla $T_{n}\leq t<T_{{n+1}}$ . Następujący algorytm formalizuje opisany powyżej sposób generacji.

Gen ${\hat{X}}_{0}\sim q(\cdot)$ ;

$T_{0}:=0$ ;

for $i:=1$ to $\infty$ do

begin

Gen $W_{{i}}\sim{\rm Ex}(Q({\hat{X}}_{{i-1}}))$ ;

$T_{{i}}:=T_{{i-1}}+W_{{i}}$ ;

Gen ${\hat{X}}_{{i}}\sim{\hat{P}}({\hat{X}}_{{i-1}},\;\cdot\;)$ ; { ${\hat{P}}$ dane wzorem powyżej }

end

Rzecz jasna, w praktyce trzeba ustalić sobie skończony horyzont czasowy $h$ i przerwać symulacje gdy $T_{i}>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 ${\cal X}=\{ 0,1,2,\ldots\}$ .

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\in{\cal X}$ symulujemy niezależnie jednorodny proces Poissona $\mathbb{R}^{{x}}$ o punktach skoku $R_{{1}}^{{x}}<\cdots<R_{{k}}^{{x}}<\cdots$ przy czym ten proces ma intensywność $Q(x)$ . Jeśli teraz, w drugiej fazie ${\hat{X}}_{{i-1}}=X(T_{{i-1}})=x$ , to za następny moment skoku wybierzemy najbliższy punkt procesu $\mathbb{R}^{{x}}$ na prawo od $T_{{i-1}}$ , czyli $T_{{i}}=\min\{ R_{k}^{{x}}:R_{k}^{{x}}>T_{{i-1}}\}$ . Z własności procesu Poissona (patrz Stwierdzenie 6.2 i jego dowód) wynika, że zmienna $T_{{i}}-T_{{i-1}}$ ma rozkład wykładniczy z parametrem $Q(x)$ i metoda jest poprawna. Naprawdę nie ma nawet potrzeby generowania wszystkich procesów Poissona $\mathbb{R}^{{x}}$ . 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 $\mathbb{R}^{{*}}$ będzie procesem Poissona o intensywności $Q^{*}\geqslant\max _{x}Q(x)$ i oznaczmy jego punkty skoków przez $\{ R_{{1}}<\cdots<R_{{k}}<\cdots\}$ . Jeśli w drugiej fazie algorytmu mamy $X(R_{{i-1}})=x$ w pewnym momencie $R_{{i-1}}\in\mathbb{R}^{{*}}$ , to zaczynamy przerzedzać proces $\mathbb{R}^{{*}}$ w taki sposób, aby otrzymać proces o intensywności $Q(x)\leqslant Q^{*}$ . Dla każdego punktu $R_{j}$ , $j\geqslant i$ powinniśmy rzucić monetą i z prawdopodobieństwem $(1-Q(x))/Q^{*}$ ten punkt usunąć. Ale jeśli usuniemy punkt $R_{i}$ to znaczy, że w tym punkcie nie ma skoku, czyli $X(R_{{i}})=x$ . Jeśli punkt $R_{i}$ zostawimy, to wykonujemy skok, losując kolejny stan z prawdopodobieństwem ${\hat{P}}(x,\cdot)$ . W rezultacie, stan $X(R_{{i}})$ jest wylosowany z rozkładu prawdopodobieństwa

$\qquad{\tilde{P}}(x,y)=\begin{cases}Q(x,y)/Q^{*}&\text{jeśli $x\not=y$;}\\ 1-Q(x)/Q^{*}&\text{jeśli $x=y$}.\end{cases}$

(7.7)

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

Odpowiada temu następujący dwufazowy algorytm.

Faza I

Gen $\mathbb{R}\sim{\rm Poiss}(Q^{*})$ { proces Poissona }

Faza II

Gen ${\tilde{X}}_{0}\sim q(\cdot)$ ;

for $i:=1$ to $\infty$ do Gen ${\tilde{X}}_{{i}}\sim{\tilde{P}}({\hat{X}}_{{i-1}},\;\cdot\;)$ ; { ${\tilde{P}}$ dane wzorem powyżej }

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 ${\cal X}=\{ 0,1,2,\ldots\}^{d}$ ; stanem układu jest układ liczb $x=(x(1),\ldots,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 $\max _{x}Q(x)=\infty$ , 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.