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=P⁢igła przetnie którąś z prostych=2⁢lπ⁢d.

Buffon zauważył, że tę prostą obserwację można wykorzystać do…obliczania liczby π 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

⁢p⌃n⁢=liczba sukcesówliczba doświadczeń⁢⁢=1n⁢∑i=1nI⁢w i-tym doświadczeniu igła przecięła prostą⁢⁢

Wielkość 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 π przyjąć

πn=2⁢lp⌃n⁢d.

Statystyk powiedziałby, że zmienna losowa πn jest estymatorem liczby π. 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 V,E będzie grafem skierowanym spójnym. Krawędzie reprezentują ,,połączenia” pomiędzy wierzchołkami. Krawędź e∈E, niezależnie od pozostałych, ulega awarii ze znanym prawdopodobieństwem pe. W rezultacie powstaje losowy podzbiór C⊆E sprawnych połączeń o rozkładzie prawdopodobieństwa

P⁢C=∏e∉Cpe⁢∏e∈C1-pe.

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

e∉C jeśli ⁢Ue≤pe;e∈C jeśli ⁢Ue>pe.

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

θ=P⁢ w zbiorze C istnieje ścieżka z v0 do v1.

Generujemy niezależnie n kopii C1,…,Cn zbioru C. Nieznane prawdopodobieństwo θ przybliżamy przez odpowiednik próbkowy:

θ⌃n=1n⁢∑i=1nI⁢ w zbiorze Ci istnieje ścieżka z v0 do v1.

Przykład 1.3 (Skomplikowana całka)

Załóżmy, że X1,…,Xm są niezależnymi zmiennymi losowymi o jednakowym rozkładzie N⁢0,1. Niech

R=maxk=1m⁡∑i=1kXi-mink=1m⁡∑i=1kXi.

Chcemy obliczyć dystrybuantę zmiennej losowej R,

H⁢x=PR≤x.

Zauważmy, że z definicji,

H⁢x=∫⋯⁢∫Ωm2⁢π-m/2⁢exp⁡-12⁢∑i=1mxi2⁢d⁢x1⁢⋯⁢d⁢xm,

gdzie Ωm=x1,…,xm:maxk=1m⁡∑i=1kxi-mink=1m⁡∑i=1kxi. W zasadzie, jest to więc zadanie obliczenia całki. Jednak skomplikowany kształt wielowymiarowego zbioru Ω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 R1,…,Rn o rozkładzie takim samym jak R. Wystarczy teraz policzyć ile spośród tych zmiennych jest ≤x:

H⌃n⁢x=1n⁢∑i=1nIRi≤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 D będzie podzbiorem kraty całkowitoliczbowej Z2. Oznaczmy przez ∂⁡D brzeg tego zbioru, zdefiniowany następująco:

⁢x,y∈∂⁡D⁢≡x,y∉D⁢ i x+1,y∈D⁢ lub ⁢x-1,y∈D⁢⁢⁢⁢⁢⁢⁢ lub ⁢x,y+1∈D⁢ lub ⁢x,y-1∈D.

Powiemy, że u:D∪∂⁡D→R jest funkcją harmoniczną, jeśli dla każdego punktu x,y∈D,

u⁢x,y=14⁢∑x′,y′∈∂⁡x,yu⁢x′,y′,

gdzie sumowanie rozciąga się na 4 punkty x′,y′ sąsiadujące z x,y, to znaczy ∂⁡x,y=x+1,y,x-1,y,x,y+1,x,y-1.

Mamy daną funkcję na brzegu: u¯:∂⁡D→R. Zadanie polega na skonstruowaniu jej rozszerzenia harmonicznego, to znaczy takiej funkcji harmonicznej u:D∪∂⁡D→R, że u⁢x,y=u¯⁢x,y dla x,y∈∂⁡D.

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

⁢⁢X0,Y0=x,y;⁢⁢Xk+1,Yk+1=Xk,Yk+ξk+1,ηk+1,

gdzie ξk,ηk są niezależnymi wektorami losowymi o jednakowym rozkładzie prawdopodobieństwa:

P((ξk,ηk)=(0,1))=P((ξk,ηk)=(0,-1))=P((ξk,ηk)=(1,0))=P((ξk,ηk)=(-1,0))=14.

Błądzimy tak długo, aż natrafimy na brzeg obszaru. Formalnie, określamy moment zatrzymania T=min⁡k:Xk,Yk∈∂⁡D. Łatwo zauważyć, że funkcja

u(x,y):=E[u¯(XT,YT)|(X0,Y0)=(x,y)]

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

u(x,y)=14∑x′,y′∈∂⁡x,yE[u¯(XT,YT)|(X1,Y1)=(x′,y′)].

Wystarczy teraz spostrzeżenie, że rozkład zmiennej losowej u⁢XT,YT pod warunkiem X1,Y1=x′,y′ jest taki sam jak pod warunkiem X0,Y0=x′,y′, 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 u¯⁢XT,YT”

• 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ń:

Przykład 1.5 (Problem ,,plecakowy”)

Załóżmy, że a=a1,…,am⊤ jest wektorem o współrzędnych naturalnych (ai∈1,2,…) i b∈1,2,…. Rozważamy wektory x=x1,…,xm⊤ o współrzędnych zero-jedynkowych (xi∈0,1). Interesuje nas liczba rozwiązań nierówności

x⊤⁢a=∑i=1mxi⁢ai≤b,

a więc liczność X⁢b zbioru X⁢b=x∈0,1m:x⊤⁢a≤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,1m, to PX∈X⁢b=X⁢b/2m. 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,1m 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,…,1⊤ i b=m/3, to prawdopodobieństwo jest ≤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 ≥em/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 a1≤⋯≤am. Niech b0=0 oraz bj=min⁡b,∑i=1jai. Rozważmy ciąg zadań plecakowych ze zmniejszającą się prawą stroną nierówności, równą kolejno b=bm,bm-1,…,b1,b0=0. Niech więc X⁢bj=x∈0,1m:x⊤⁢a≤bj. Zachodzi następujący ,,wzór teleskopowy”:

X⁢b=X⁢bm=X⁢bmX⁢bm-1⋅X⁢bm-1X⁢bm-2⁢⋯⁢X⁢b1X⁢b0⋅X⁢b0.

Oczywiście, |X(b0)=1|, a więc możemy uznać, że zadanie sprowadza się do obliczenia ilorazów X⁢bj-1/X⁢bj (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⁢bj, to moglibyśmy postępować w dobrze już znany sposób: liczyć ,,sukcesy” polegające na wpadnięciu w zbiór X⁢bj-1⊆X⁢bj. Rzecz jasna, możemy losować z rozkładu jednostajnego na kostce 0,1m i eliminować punkty poza zbiorem X⁢bj, 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 ≥em/18 losowań?

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

Zaczynamy błądzenie w punkcie 0=0,…,0⊤. 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 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:

PX⁢n=x→1X⁢b⁢n→∞