Zagadnienia

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.

Treść automatycznie generowana z plików źródłowych LaTeXa za pomocą oprogramowania wykorzystującego LaTeXML.

Projekt współfinansowany przez Unię Europejską w ramach Europejskiego Funduszu Społecznego.

Projekt współfinansowany przez Ministerstwo Nauki i Szkolnictwa Wyższego i przez Uniwersytet Warszawski.