Zagadnienia

1. Wprowadzenie

1.1. Od autora

Jest kilka ważnych powodów, dla których warto się zająć symulacjami stochastycznymi:

  • Symulacje stochastyczne są prostym sposobem badania zjawisk losowych.

  • Ściśle związane z symulacjami stochastycznymi są metody obliczeniowe nazywane ,,Monte Carlo” (MC). Polegają one na wykorzystaniu ,,sztucznie generowanej” losowości w celu rozwiązania zadań deterministycznych. Metody MC są proste i skuteczne. Dla pewnych problemów MC jest jedynym dostępnym narzędziem obliczeniowym. Dla innych problemów MC jest co prawda mniej efektywne od metod numerycznych, ale za to dużo łatwiejsze!

  • W moim przekonaniu symulacje stochastyczne są wspaniałą pomocą przy nauce rachunku prawdopodobieństwa. Pozwalają lepiej ,,zrozumieć losowość”.

  • Symulacje stochastyczne są dostępne dla każdego. W szczególności, ,,otoczenie” R, które stanowi naprawdę potężne narzędzie, jest rozpowszechniane za darmo!

Jest wreszcie powód najważniejszy:

  • Symulacje stochastyczne są świetną zabawą!

Literatura na temat symulacji stochastycznych jest bardzo obszerna. Godna polecenia jest książka Zielińskiego i Wieczorkowskiego [23], poświęcona w całości generatorom zmiennych losowych. Przedstawia ona bardziej szczegółowo zagadnienia, odpowiadające Rozdziałom 2–4 niniejszego skryptu i zawiera materiał, który zdecydowałem się pominąć: wytwarzanie ,,liczb losowych” o rozkładzie jednostajnym i testowanie generatorów. Podobne zagadnienia są przedstawione trochę w innym stylu w monografii Ripleya [18], która również zawiera wstęp do metod Monte Carlo. Zaawansowane wykłady można znależć w nowoczesnych monografiach Asmussena i Glynna [2], Liu [15], Roberta i Caselli [19]. Pierwsza z nich jest zorientowana bardziej na wyniki teoretyczne, zaś druga bardziej na zastosowania. Świetnym wstępem do metod MCMC są prace Geyera [7] i [8]. Teoria łańcuchów Markowa z uwzględnieniem zagadnień istotnych dla MCMC jest przystępnie przedstawiona w książce Brémaud [4]. Podstawy teoretycznej analizy zrandomizowanych algorytmów (tematy poruszane w Rozdziale 15 skryptu) są znakomicie przedstawione w pracach Jerruma i Sinclaira [11] oraz Jerruma [12].

1.2. Przykłady

Zacznę od kilku przykładów zadań obliczeniowych, które można rozwiązywać symulując losowość. Wybrałem przykłady najprostsze, do zrozumienia których wystarcza zdrowy rozsądek i nie potrzeba wielkiej wiedzy.

Przykład 1.1 (Igła Buffona, 1777)

Podłoga jest nieskończoną płaszczyzną, podzieloną równoległymi prostymi na ,,deski” szerokości d. Rzucamy ,,losowo” igłę o długości l. Elementarne rozważania prowadzą do wniosku, że

p=\mathbb{P}(\text{igła przetnie którąś z prostych})=\frac{2l}{\pi d}.

Buffon zauważył, że tę prostą obserwację można wykorzystać do…obliczania liczby \pi metodą ,,statystyczną”. Powtórzmy nasze doświadczenie niezależnie n razy. Oczywiście, mamy do czynienia ze schematem Bernoulliego, w którym ,,sukcesem” jest przecięcie prostej. Niech

\begin{split}{\hat{p}_{n}}&=\frac{\text{liczba sukcesów}}{\text{liczba doświadczeń}}\\
&=\frac{1}{n}\sum _{{i=1}}^{n}\mathbb{I}(\text{w $i$-tym doświadczeniu igła przecięła prostą})\\
\end{split}

Wielkość {\hat{p}_{n}} jest empirycznym odpowiednikiem prawdopodobieństwa p i powinna przybliżać to prawdopodobieństwo, przynajmniej dla dużych n. W takim razie, możemy za przybliżenie liczby \pi przyjąć

{\pi _{n}}=\frac{2l}{{\hat{p}_{n}}d}.

Statystyk powiedziałby, że zmienna losowa {\pi _{n}} jest estymatorem liczby \pi. To wszystko jest bardzo proste, ale parę kwestii wymaga uściślenia. Jak duża ma być liczba powtórzeń, żeby przybliżenie było odpowiednio dokładne? Ale przecież mamy do czynienia ze ,,ślepym losem”! Czyż pech nie może sprawić, że mimo dużej liczby doświadczeń przybliżenie jest kiepskie? Czy odpowiednio dobierając l i d możemy poprawić dokładność? A może da się zaprojektować lepsze doświadczenie?

Przykład 1.2 (Sieć zawodnych połączeń)

Niech {\cal V},{\cal E} będzie grafem skierowanym spójnym. Krawędzie reprezentują ,,połączenia” pomiędzy wierzchołkami. Krawędź e\in{\cal E}, niezależnie od pozostałych, ulega awarii ze znanym prawdopodobieństwem p_{e}. W rezultacie powstaje losowy podzbiór C\subseteq{\cal E} sprawnych połączeń o rozkładzie prawdopodobieństwa

P(C)=\prod _{{e\not\in C}}p_{e}\prod _{{e\in C}}(1-p_{e}).

Jest jasne, jak można symulować to zjawisko: dla każdej krawędzi e ,,losujemy” zmienną U_{e} o rozkładzie równomiernym na [0,1] i przyjmujemy

\begin{cases}e\not\in C&\text{ jeśli }U_{e}\le p_{e};\\
e\in C&\text{ jeśli }U_{e}>p_{e}.\\
\end{cases}

Powiedzmy, że interesuje nas możliwość znalezienia ścieżki (ciągu krawędzi) wiodącej z ustalonego wierzchołka v_{0} do innego wierzchołka v_{1}. Niech

\theta=\mathbb{P}(\text{ w zbiorze $C$ istnieje ścieżka z $v_{0}$ do $v_{1}$}).

Generujemy niezależnie n kopii C_{1},\ldots,C_{n} zbioru C. Nieznane prawdopodobieństwo \theta przybliżamy przez odpowiednik próbkowy:

{\hat{\theta}_{n}}=\frac{1}{n}\sum _{{i=1}}^{n}\mathbb{I}(\text{ w zbiorze $C_{i}$ istnieje ścieżka z $v_{0}$ do $v_{1}$}).
Przykład 1.3 (Skomplikowana całka)

Załóżmy, że X_{1},\ldots,X_{m} są niezależnymi zmiennymi losowymi o jednakowym rozkładzie {\rm N}(0,1). Niech

R=\max _{{k=1}}^{m}\sum _{{i=1}}^{k}X_{i}-\min _{{k=1}}^{m}\sum _{{i=1}}^{k}X_{i}.

Chcemy obliczyć dystrybuantę zmiennej losowej R,

H(x)=\mathbb{P}(R\leq x).

Zauważmy, że z definicji,

H(x)={\int\cdots\int}_{{\Omega _{m}}}(2\pi)^{{-m/2}}\exp\left[-\frac{1}{2}\sum _{{i=1}}^{m}x_{i}^{2}\right]{\rm d}x_{1}\cdots{\rm d}x_{m},

gdzie \Omega _{m}=\left\{(x_{1},\ldots,x_{m}):\max _{{k=1}}^{m}\sum _{{i=1}}^{k}x_{i}-\min _{{k=1}}^{m}\sum _{{i=1}}^{k}x_{i}\right\}. W zasadzie, jest to więc zadanie obliczenia całki. Jednak skomplikowany kształt wielowymiarowego zbioru \Omega _{m} powoduje, że zastosowanie standardowych metod numerycznych jest utrudnione. Z kolei symulowanie zmiennej losowej R jest bardzo łatwe, wprost z definicji. Można wygenerować wiele niezależnych zmiennych R_{1},\ldots,R_{n} o rozkładzie takim samym jak R. Wystarczy teraz policzyć ile spośród tych zmiennych jest \leq x:

\hat{H}_{n}(x)=\frac{1}{n}\sum _{{i=1}}^{n}\mathbb{I}(R_{i}\leq x).

Schemat symulacyjnego obliczania prawdopodobieństwa jest taki sam jak w dwu poprzednich przykładach. Podobnie zresztą jak w tamtych przykładach, podstawowy algorytm można znacznie ulepszać, ale dyskusję na ten temat odłóżmy na póżniej.

Przykład 1.4 (Funkcja harmoniczna)

Następujące zadanie jest dyskretnym odpowiednikiem sławnego zagadnienia Dirichleta. Niech {\cal D} będzie podzbiorem kraty całkowitoliczbowej \mathbb{Z}^{2}. Oznaczmy przez \partial{\cal D} brzeg tego zbioru, zdefiniowany następująco:

\begin{split}(x,y)\in\partial{\cal D}&\equiv(x,y)\not\in{\cal D}\text{ i }\big[(x+1,y)\in{\cal D}\text{ lub }(x-1,y)\in{\cal D}\\
&\qquad\qquad\qquad\quad\text{ lub }(x,y+1)\in{\cal D}\text{ lub }(x,y-1)\in{\cal D}\big].\end{split}

Powiemy, że u:{\cal D}\cup\partial{\cal D}\to\mathbb{R} jest funkcją harmoniczną, jeśli dla każdego punktu (x,y)\in{\cal D},

u(x,y)=\frac{1}{4}\sum _{{(x^{\prime},y^{\prime})\in\partial(x,y)}}u(x^{\prime},y^{\prime}),

gdzie sumowanie rozciąga się na 4 punkty (x^{\prime},y^{\prime}) sąsiadujące z (x,y), to znaczy \partial(x,y)=\{(x+1,y),(x-1,y),(x,y+1),(x,y-1)\}.

Mamy daną funkcję na brzegu: \bar{u}:\partial{\cal D}\to\mathbb{R}. Zadanie polega na skonstruowaniu jej rozszerzenia harmonicznego, to znaczy takiej funkcji harmonicznej u:{\cal D}\cup\partial{\cal D}\to\mathbb{R}, że u(x,y)=\bar{u}(x,y) dla (x,y)\in\partial{\cal D}.

Wyobraźmy sobie błądzenie losowe po kracie, startujące w punkcie (x,y)\in{\cal D}. Formalnie jest to ciąg losowych punktów określonych następująco:

\begin{split}&(X_{0},Y_{0})=(x,y);\\
&(X_{{k+1}},Y_{{k+1}})=(X_{k},Y_{k})+(\xi _{{k+1}},\eta _{{k+1}}),\end{split}

gdzie (\xi _{k},\eta _{k}) są niezależnymi wektorami losowymi o jednakowym rozkładzie prawdopodobieństwa:

\begin{split}\mathbb{P}\left((\xi _{k},\eta _{k})=(0,1)\right)&=\mathbb{P}\left((\xi _{k},\eta _{k})=(0,-1)\right)\\
=\mathbb{P}\left((\xi _{k},\eta _{k})=(1,0)\right)&=\mathbb{P}\left((\xi _{k},\eta _{k})=(-1,0)\right)=\frac{1}{4}.\end{split}

Błądzimy tak długo, aż natrafimy na brzeg obszaru. Formalnie, określamy moment zatrzymania T=\min\{ k:(X_{k},Y_{k})\in\partial{\cal D}\}. Łatwo zauważyć, że funkcja

u(x,y):=\mathbb{E}[\bar{u}(X_{T},Y_{T})|(X_{0},Y_{0})=(x,y)]

jest rozwiązaniem zagadnienia! Istotnie, ze wzoru na prawdopodobieństwo całkowite wynika, że

u(x,y)=\frac{1}{4}\sum _{{(x^{\prime},y^{\prime})\in\partial(x,y)}}\mathbb{E}[\bar{u}(X_{T},Y_{T})|(X_{1},Y_{1})=(x^{\prime},y^{\prime})].

Wystarczy teraz spostrzeżenie, że rozkład zmiennej losowej u(X_{T},Y_{T}) pod warunkiem (X_{1},Y_{1})=(x^{\prime},y^{\prime}) jest taki sam jak pod warunkiem (X_{0},Y_{0})=(x^{\prime},y^{\prime}), bo błądzenie ,,rozpoczyna się na nowo”.

Algorytm Monte Carlo obliczania u(x,y) oparty na powyższym spostrzeżeniu został wynaleziony przez von Neumanna i wygląda następująco:

  • Powtórz wielokrotnie, powiedzmy n razy, niezależnie doświadczenie:

  • ,,błądź startując startując z (x,y) aż do brzegu; oblicz \bar{u}(X_{T},Y_{T})

  • Uśrednij wyniki n doświadczeń.

Dla bardziej formalnego zapisu algorytmu będę się posługiwał pseudo-kodem, który wydaje się zrozumiały bez dodatkowych objaśnień:

{ 'Gen' oznacza 'Generuj' }
U:=0;
for j=1 to n
  begin
  (X,Y):=(x,y);
  while (X,Y)\in{\cal D}
    begin
    Gen (\xi,\eta)\sim{\rm U}\{(0,1),(0,-1),(1,0),(-1,0)\}; [ rozkład jednostajny na zbiorze 4-punktowym ]
    (X,Y):=(X,Y)+(\xi,\eta)
    end
    U:=U+\bar{u}(X,Y)
  end
U:=U/n
Przykład 1.5 (Problem ,,plecakowy”)

Załóżmy, że a=(a_{1},\ldots,a_{m})^{\top} jest wektorem o współrzędnych naturalnych (a_{i}\in\{ 1,2,\ldots\}) i b\in\{ 1,2,\ldots\}. Rozważamy wektory x=(x_{1},\ldots,x_{m})^{\top} o współrzędnych zero-jedynkowych (x_{i}\in\{ 0,1\}). Interesuje nas liczba rozwiązań nierówności

x^{\top}a=\sum _{{i=1}}^{m}x_{i}a_{i}\leq b,

a więc liczność |{\cal X}(b)| zbioru {\cal X}(b)=\{ x\in\{ 0,1\}^{m}:x^{\top}a\leq b\}. Czytelnik domyśla się z pewnością, skąd nazwa problemu. Dokładne obliczenie liczby rozwiązań jest trudne. Można spróbować zastosować prostą metodę Monte Carlo. Jeśli zmienna losowa X ma rozkład jednostajny na przestrzeni \{ 0,1\}^{m}, to \mathbb{P}(X\in{\cal X}(b))=|{\cal X}(b)|/2^{m}. Wystarczy zatem oszacować to prawdopodobieństwo tak jak w trzech pierwszych przykładach w tym rozdziale. Generowanie zmiennej losowej o rozkłdzie jednostajnym na przestrzeni \{ 0,1\}^{m} sprowadza się do przeprowadzenia m rzutów monetą i jest dziecinnie proste. Na czym więc polega problem? Otóż szacowane prawdopodobieństwo może być astronomicznie małe. Dla, powiedzmy a=(1,\ldots,1)^{\top} i b=m/3, to prawdopodobieństwo jest \leq{\rm e}^{{-m/18}} (proszę się zastanowić jak uzasadnić tę nierówność). Przeprowadzanie ciągu doświadczeń Bernoulliego z takim prawdopodobieństwem sukcesu jest bardzo nieefektywne – na pierwszy sukces oczekiwać będziemy średnio \geq{\rm e}^{{m/18}}, co dla dużych m jest po prostu katastrofalne.

Metoda, którą naszkicuję należy do rodziny algorytmów MCMC (Monte Carlo opartych na łańcuchach Markowa). Algorytm jest raczej skomplikowany, ale o ile mi wiadomo jest najefektywniejszym ze znanych sposobów rozwiązania zadania. Bez straty ogólności możemy przyjąć, że a_{1}\leq\cdots\leq a_{m}. Niech b_{0}=0 oraz b_{j}=\min\{ b,\sum _{{i=1}}^{j}a_{i}\}. Rozważmy ciąg zadań plecakowych ze zmniejszającą się prawą stroną nierówności, równą kolejno b=b_{m},b_{{m-1}},\ldots,b_{1},b_{0}=0. Niech więc {\cal X}(b_{j})=\{ x\in\{ 0,1\}^{m}:x^{\top}a\leq b_{j}\}. Zachodzi następujący ,,wzór teleskopowy”:

|{\cal X}(b)|=|{\cal X}(b_{m})|=\frac{|{\cal X}(b_{{m}})|}{|{\cal X}(b_{{m-1}})|}\cdot\frac{|{\cal X}(b_{{m-1}})|}{|{\cal X}(b_{{m-2}})|}\cdots\frac{|{\cal X}(b_{{1}})|}{|{\cal X}(b_{{0}})|}\cdot|{\cal X}(b_{{0}})|.

Oczywiście, |{\cal X}(b_{{0}})=1|, a więc możemy uznać, że zadanie sprowadza się do obliczenia ilorazów {|{\cal X}(b_{{j-1}})|}/{|{\cal X}(b_{{j}})|} (pomijamy w tym miejscu subtelności związane z dokładnością obliczeń). Gdybyśmy umieli efektywnie generować zmienne losowe o rozkładzie jednostajnym na przestrzeni X(b_{{j}}), to moglibyśmy postępować w dobrze już znany sposób: liczyć ,,sukcesy” polegające na wpadnięciu w zbiór X(b_{{j-1}})\subseteq X(b_{{j}}). Rzecz jasna, możemy losować z rozkładu jednostajnego na kostce \{ 0,1\}^{m} i eliminować punkty poza zbiorem X(b_{{j}}), ale w ten sposób wpadlibyśmy w pułapkę, od której właśnie uciekamy: co zrobić jeśli przez sito eliminacji przechodzi jedno na \geq{\rm e}^{{m/18}} losowań?

Opiszemy pewne wyjście, polegające na zastosowaniu błądzenia losowego po zbiorze X(b_{{j}}). Dla ustalenia uwagi przyjmijmy j=m, czyli b_{j}=b. Losowy ciąg punktów X(0),X(1),\ldots,X(n),\ldots generujemy rekurencyjnie w taki sposób:

X(0):=0;
for n=0 to \infty
  begin
  X:=X(n); [ gdzie X=(X_{1},\ldots,X_{m})
  Gen I\sim{\rm U}\{ 1,\ldots,m\}; [ losujemy indeks do zmiany ]
  Y:=XY_{I}=1-X_{I};
  if Y^{\top}a\leq b then X(n+1):=Y
       else X(n+1):=X
  end
end

Zaczynamy błądzenie w punkcie 0=(0,\ldots,0)^{\top}. Jeśli w chwili n jesteśmy w punkcie X, to próbujemy przejść do nowego punktu Y, utworzonego przez zmianę jednej, losowo wybranej współrzędnej punktu X (0 zamieniamy na 1 lub z 1 na 0; pozostałe współrzędne zostawiamy bez zmian). Jeśli ,,proponowany punkt Y nie wypadł” z rozważanej przestrzeni {\cal X}(b), to przechodzimy do punktu Y. W przeciwnym wypadku stoimy w punkcie X.

Rzecz jasna, generowane w ten sposób zmienne losowe X(n) nie mają dokładnie rozkładu jednostajnego ani tym bardziej nie są niezależne. Jednak dla dużych n zmienna X(n) ma w przybliżeniu rozkład jednostajny:

\mathbb{P}(X(n)=x)\to\frac{1}{|{\cal X}(b)|}\qquad(n\to\infty)

dla każdego x\in{\cal X}(b). Co więcej, z prawdopodobieństwem 1 jest prawdą, że

\frac{1}{n}\sum _{{i=1}}^{n}\mathbb{I}(X(i)=x)\to\frac{1}{|{\cal X}(b)|}\qquad(n\to\infty).

Innymi słowy, ciąg X(n) spełnia prawo wielkich liczb – i może być użyty do szacowania liczności podzbiorów przestrzeni {\cal X}(b) zamiast trudnego do symulowania ciągu niezależnych zmiennych o rozkładzie jednostajnym.

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.