Zagadnienia

8.1 Losowanie istotne
8.2 Efektywność estymatorów MC
8.3 Ważona eliminacja

8. Algorytmy Monte Carlo I. Obliczanie całek

Duża część zastosowań metod MC wiąże się z obliczaniem całek lub sum. Typowe zadanie polega na obliczeniu wartości oczekiwanej

$\qquad\theta=\mathbb{E}_{{\pi}}f(X)=\int _{{{\cal X}}}f(x)\pi({\rm d}x),$

gdzie $X$ jest zmienną losową o rozkładzie prawdopodobieństwa $\pi$ na przestrzeni ${\cal X}$ , zaś $f:{\cal X}\to\mathbb{R}$ . Zazwyczaj ${\cal X}$ jest podzbiorem wielowymiarowej przestrzeni euklidesowej lub jest zbiorem skończonym ale bardzo licznym (na przykład zbiorem pewnych obiektów kombinatorycznych). Jeśli ${\cal X}\subseteq\mathbb{R}^{{d}}$ i rozkład $\pi$ jest opisany przez gęstość $p$ to całka określająca wartość oczekiwaną jest zwykłą całką Lebesgue'a. W tym przypadku zatem możemy napisać

$\qquad\theta=\int _{{{\cal X}}}f(x)p(x){\rm d}x.$

Równie ważny w zastosowaniach jest przypadek dyskretnej przestrzeni ${\cal X}$ . Jeśli $p(x)=\mathbb{P}(X=x)$ , to

$\qquad\theta=\sum _{{x\in{\cal X}}}f(x)p(x).$

W dalszym ciągu tego rozdziału utożsamiamy rozkład $\pi$ z funkcją $p$ . Zwróćmy uwagę, że przedstawiamy tu $\theta$ jako wartość oczekiwaną tylko dla wygody oznaczeń; w istocie każda całka, suma, (a także prawdopodobieństwo zdarzenia losowego) jest wartością oczekiwaną.

Na pierwszy rzut oka nie widać czym polega problem! Sumowanie wykonuje każdy kalkulator, całki się sprawnie oblicza numerycznie. Ale nie zawsze. Metody MC przydają się w sytuacjach gdy spotyka się z następującymi trudnościami (lub przynajmniej którąś z nich).

Przestrzeń ${\cal X}$ jest ogromna. To znaczy wymiar $d$ jest bardzo duży lub skończona przestrzeń zawiera astronomicznie dużą liczbę punktów.
Rozkład prawdopodobieństwa $\pi$ jest ,,skupiony” w małej części ogromnej przestrzeni ${\cal X}$ .
Nie ma podstaw do zakładania, że funkcja $f$ jest w jakimkolwiek sensie ,,gładka” (co jest warunkiem stosowania standardowych numerycznych metod całkowania).
Gęstość $p$ rozkładu $\pi$ jest znana tylko z dokładnością do stałej normującej. Innymi słowy, umiemy obliczać ${p}^{\prime}(x)=Z{p}(x)$ ale nie znamy stałej $Z=\int{p}^{\prime}(x){\rm d}x$ . Czsem zadanie polega właśnie na obliczeni tej stałej ( $Z$ jest nazwane ,,funkcją podziału” lub sumą statystyczną).

Prosta metoda MC (po angielsku nazywana bardziej brutalnie: Crude Monte Carlo, czyli CMC) nasuwa się sama. Należy wygenerować $n$ niezależnych zmiennych losowych $X_{{1}},\ldots,X_{{n}}$ o jednakowym rozkładzie $\pi$ i za estymator wartości oczekiwanej wziąć średnią z próbki,

$\hat{\theta}_{n}=\hat{\theta}_{n}^{{\rm CMC}}=\sum _{{i=1}}^{n}f(X_{{i}}).$

Mocne Prawo Wielkich Liczb gwarantuje, że $\hat{\theta}_{n}\to\theta$ prawie na pewno, gdy $n\to\infty$ . W terminologii statystycznej, $\hat{\theta}_{n}$ jest mocno zgodnym estymatorem obliczanej wielkości. Zadanie wydaje się rozwiązane. Są jednak dwa zasadnicze kłopoty.

Estymator $\hat{\theta}_{n}$ skonstruowany metodą CMC może się zbliżać do $\theta$ przerażająco wolno.
Co robić, jeśli nie umiemy losować z rozkładu $\pi$ ?

8.1. Losowanie istotne

Zadziwiająco skutecznym sposobem na oba przedstawione wyżej kłopoty jest losowanie istotne (Importance Sampling, w skrócie IS). Przypuśćmy, że umiemy losować z rozkładu o gęstości $q$ . Zauważmy, że

$\qquad\theta=\int _{{{\cal X}}}\frac{p(x)}{q(x)}f(x)q(x){\rm d}x=\mathbb{E}_{{q}}\frac{p(X)}{q(X)}f(X)=\mathbb{E}_{{q}}w(X)f(X),$

gdzie $w(x)={p}(x)/q(x)$ . Piszemy tu wzory dla całek ale oczywiście dla sum jest tak samo. Niech $X_{{1}},\ldots,X_{{n}}$ będą niezależnymi zmiennymi losowymi o jednakowym rozkładzie $g$ ,

$\qquad\hat{\theta}_{n}=\hat{\theta}_{n}^{{\rm IS1}}=\sum _{{i=1}}^{n}W_{{i}}f(X_{{i}}),$

(8.1)

gdzie

$\qquad W_{{i}}=w(X_{{i}})=\frac{{p}(X_{{i}})}{q(X_{{i}})}$

traktujemy jako wagi wylosowanych punktów $X_{{i}}$ . Z tego co wyżej powiedzieliśmy wynika, że $\mathbb{E}_{{q}}\hat{\theta}_{n}=\theta$ oraz $\hat{\theta}_{n}\to\theta$ prawie na pewno, gdy $n\to\infty$ . Mamy więc estymator nieobciążony i zgodny. Milcząco założyliśmy, że $q(X_{{i}})>0$ w każdym wylosowanym punkcie, czyli, że $\{ x:{p}(x)>0\}\subseteq\{ x:q(x)>0\}$ . Z tym zwykle nie ma wielkiego kłopotu. Ponadto musimy założyć, że umiemy obliczać funkcję $w$ w każdym wylosowanym punkcie. Jeśli znamy tylko ${p}^{\prime}(x)=Z{p}(x)$ a nie znamy stałej $Z$ to jest kłopot. Możemy tylko obliczyć $W^{\prime}_{i}={p}^{\prime}(X_{i})/q(X_{{i}})=ZW_{{i}}$ . Używamy zatem nieco innej postaci estymatora IS, mianowicie

$\qquad\hat{\theta}_{n}=\hat{\theta}_{n}^{{\rm IS2}}=\frac{\sum _{{i=1}}^{n}W_{{i}}f(X_{{i}})}{\sum _{{i=1}}^{n}W_{{i}}}.$

(8.2)

Ten estymator wygląda dziwnie, bo w mianowniku występuje estymator jedynki ale chodzi o to, że we wzorze (8.2) w przeciwieństwie do (8.1) możemy zastąpić $W_{{i}}$ przez $W^{\prime}_{{i}}$ , bo nieznany czynnik $Z$ skraca się. Możemy jeszcze zapisać nasz estymator w zgrabnej formie

$\qquad\hat{\theta}_{n}^{{\rm IS2}}=\sum _{{i=1}}^{n}\tilde{W}_{{i}}f(X_{{i}}),$

gdzie $\tilde{W}_{i}=W_{{i}}/\sum _{{j}}W_{{j}}$ są unormowanymi wagami. Trzeba jednak pamiętać, że to ,,unormowanie” polega na podzieleniu wag przez zmienną losową. W rezultacie estymator ${\rm IS2}$ jest obciążony. Jest jednak mocno zgodny, bo licznik we wzorze (8.2) po podzieleniu przez $n$ zmierza do $\theta$ a mianownik do 1, prawie na pewno. Estymator ${\rm IS2}$ jest znacznie częściej używany niż ${\rm IS1}$ . Oprócz uniezależnienia się od stałej normującej, jest jeszcze inny powód. Okazuje się, że (pomimo obciążenia) estymator ${\rm IS2}$ może być bardziej efektywny. Sens tego pojęcia wyjaśnimy w następnym podrozdziale. Zarówno dobór przykładów jak i ogólny tok rozważań jest zapożyczony z monografii Liu [15].

8.2. Efektywność estymatorów MC

Naturalne jest żądanie, aby estymator był mocno zgodny, $\hat{\theta}_{n}\to\theta$ prawie na pewno przy $n\to\infty$ . Znaczy to, że zbliżymy się dowolnie blisko aproksymowanej wielkości, gdy tylko dostatecznie długo przedłużamy symulacje. Chcielibyśmy jednak wiedzieć coś konkretniejszego, oszacować jak szybko błąd aproksymacji $\hat{\theta}_{n}-\theta$ maleje do zera i jakie $n$ wystarczy do osiągnięcia wystarczającej dokładności. Przedstawimy teraz najprostsze i najczęściej stosowane w praktyce podejście, oparte na tzw. asymptotycznej wariancji i konstrukcji asymptotycznych przedziałów ufności. W typowej sytuacji estymator $\hat{\theta}_{n}$ ma następującą własność, nazywaną asymptotyczną normalnością:

$\qquad\sqrt{n}\left(\hat{\theta}_{n}-\theta\right)\to{\rm N}(0,\sigma^{2}),\qquad(n\to\infty)$

w sensie zbieżności według rozkładu. W konsekwencji, dla ustalonej liczby $a>0$ ,

$\qquad\mathbb{P}\left(|\hat{\theta}_{n}-\theta|\leqslant\frac{z\sigma}{\sqrt{n}}\right)\to\Phi(z)-\Phi(-z)=2\Phi(z)-1,\qquad(n\to\infty),$

gdzie $\Phi$ jest dystrybuantą rozkładu ${\rm N}(0,1)$ . Często nie znamy asymptotycznej wariancji $\sigma^{{2}}$ ale umiemy skonstruować jej zgodny estymator $\hat{\sigma}_{n}^{{2}}$ . Dla ustalonej ,,małej” dodatniej liczby $\alpha$ łatwo dobrać kwantyl rozkładu normalnego $z=z_{{1-\alpha/2}}$ tak, żeby $2\Phi(z)-1=1-\alpha$ . Typowo $\alpha=0.05$ i $z=z_{{0.975}}=1.9600\approx 2$ . Otrzymujemy

$\qquad\mathbb{P}\left(|\hat{\theta}_{n}-\theta|\leqslant\frac{z_{{1-\alpha/2}}\hat{\sigma}_{{n}}}{\sqrt{n}}\right)\to 1-\alpha.\quad n\to\infty.$

Ten wzór interpretuje się w praktyce tak: estymator $\hat{\theta}_{n}$ ma błąd nie przekraczający $z\hat{\sigma}_{{n}}/{\sqrt{n}}$ z prawdopodobieństwem około $1-\alpha$ . W żargonie statystycznym $1-\alpha$ jest nazywane asymptotycznym poziomem ufności.

Chociaż opisame powyżej podejście ma swoje słabe strony, to prowadzi do prostego kryterium porównywania estymatorów. Przypuśćmy, że mamy dwa estymatory $\hat{\theta}_{n}^{{\rm I}}$ i $\hat{\theta}_{n}^{{\rm II}}$ , oba asymptotycznie normalne, o asymptotycznej wariancji $\sigma _{{\rm I}}^{2}$ i $\sigma _{{\rm II}}^{6}$ odpowiednio. Przypuśćmy dalej, że dla obliczenia pierwszego z tych estymatorów generujemy $n_{{\rm I}}^{2}$ punktów, zaś dla drugiego $n_{{\rm I}}^{2}$ . Błędy obu estymatorów na tym mym poziomie istotności są ograniczone przez, odpowiednio, $z\sigma _{{\rm I}}/\sqrt{n_{{\rm I}}}$ i $z\sigma _{{\rm II}}/\sqrt{n_{{\rm II}}}$ . Przyrównując te wyrażenia do siebie dochodzimy do wniosku, że oba estymatory osiągają podobną dokładność, jeśli $n_{{\rm I}}/n_{{\rm II}}=\sigma _{{\rm II}}^{2}/\sigma _{{\rm I}}^{2}$ . Liczbę

$\qquad{\rm eff}(\hat{\theta}_{n}^{{\rm I}},\hat{\theta}_{n}^{{\rm II}})=\frac{\sigma _{{\rm II}}^{2}}{\sigma _{{\rm I}}^{2}}$

nazywamy względną efektywnością (asymptotyczną). Czasami dobrze jest wybrać za ,,naturalny punkt odniesienia” estymator CMC i zdefiniować efektywność estymatora $\hat{\theta}_{n}$ o asymptotycznej wariancji $\sigma^{{2}}$ jako

$\qquad{\rm eff}(\hat{\theta}_{n})={\rm eff}(\hat{\theta}_{n},\hat{\theta}_{n}^{{\rm CMC}})=\frac{\sigma _{{\rm CMC}}^{2}}{\sigma^{2}}.$

Mówi się też, że jeśli wygenerujemy próbkę $n$ punktów i obliczymy $\hat{\theta}_{n}$ to ,,efektywna liczność próbki” (ESS, czyli effective sample size) jest $n/{\rm eff}(\hat{\theta}_{n})$ . Tyle bowiem należałoby wygenerować punktów stosując CMC, żeby osiągnąć podobną dokładność.

Asymptotyczna normalność prostego estymatora CMC wynika wprost z Centralnego Twierdzenia Granicznego (CTG). Istotnie, jeśli generujemy niezależnie $X_{{1}},\ldots,X_{{n}}$ o jednakowym rozkładzie $\pi$ , to

$\qquad\begin{split}\sqrt{n}\left(\hat{\theta}_{n}^{{\rm CMC}}-\theta\right)&=\frac{1}{\sqrt{n}}\sum _{{i=1}}^{n}(f(X_{{i}})-\theta)\\ &\to{\rm N}(0,\sigma _{{\rm CMC}}^{2}),\qquad(n\to\infty),\end{split}$

gdzie

$\qquad\sigma _{{\rm CMC}}^{2}={\rm Var}_{{p}}f(X)=\int\left(f(x)-\theta\right)^{{2}}{p}(x){\rm d}x.$

Zupełnie podobnie, dla losowania istotnego w formie (8.1) otrzymujemy asymptotyczną normalność i przy tym

$\qquad\begin{split}\sigma _{{\rm IS1}}^{2}={\rm Var}_{{q}}w(X)f(X)&=\mathbb{E}_{q}\frac{\left(f(X){p}(X)-\theta q(X)\right)^{{2}}}{q(X)^{2}}\\ &=\int\frac{\left(f(x){p}(x)-\theta q(x)\right)^{{2}}}{q(x)}{\rm d}x.\end{split}$

Z tego wzorku widać, że estymator może mieć wariancję zero jeśli $q(x)\propto f(x){p}(x)$ . Niestety, żeby obliczyć $q(x)$ potrzebna jest znajomość współczynnika proporcjonalności, który jest równy… $\theta$ . Nie wszystko jednak stracone. Pozostaje ważna reguła heurystyczna:

Gęstość $q(x)$ należy tak dobierać, aby jej ,,kształt” był zbliżony do funkcji $f(x){p}(x)$ .

Dla drugiej wersji losowania istotnego, (8.2), asymptotyczna normalność wynika z następujących rozważań:

$\qquad\begin{split}\sqrt{n}\left(\hat{\theta}_{n}^{{\rm IS2}}-\theta\right)&=\frac{\frac{1}{\sqrt{n}}\sum _{{i=1}}^{n}(W_{{i}}f(X_{{i}})-\theta W_{{i}})}{\frac{1}{n}\sum _{{i=1}}^{n}W_{{i}}}\\ &\to{\rm N}(0,\sigma _{{\rm IS2}}^{2}),\qquad(n\to\infty),\end{split}$

gdzie

$\qquad\begin{split}\sigma _{{\rm IS2}}^{2}&={\rm Var}_{{q}}w(X)(f(X)-\theta)\\ &={\rm Var}_{{q}}w(X)f(X)-2\theta{\rm Cov}_{{q}}(w(X),w(X)f(X))+\theta^{{2}}{\rm Var}_{{q}}w(X)\\ &=\sigma _{{\rm IS1}}^{2}+\theta\left[-2{\rm Cov}_{{q}}(w(X),w(X)f(X))+\theta{\rm Var}_{{q}}w(X)\right].\end{split}$

Wyrażenie w kwadratowym nawiasie może być ujemne, jeśli jest duża dodatnia korelacja zmiennych $w(X)$ i $w(X)f(X)$ . W tej sytuacji estymator IS2 jest lepszy od IS1. Okazuje się więc rzecz na pozór paradoksalna: dzielenie przez estymator jedynki może poprawić estymator. Poza tym oba estymatory IS1 i IS2 mogą mieć mniejszą (asymptotyczną) wariancję, niż CMC. Jeśli efektywność jest większa niż 100%, to używa się czasem określenia ,,estymator superefektywny”.

Jeśli interesuje nas obliczenie wartości oczekiwanej dla wielu różnych funkcji $f$ , to warto wprowadzić uniwersalny, niezależny of $f$ wskaźnik efektywności. Niezłym takim wskaźnikiem jest ${\rm Var}_{{q}}w(X)$ . Istotnie, ,,odchylenie $\chi^{{2}}$ ” pomiędzy gęstością instrumentalną $p$ i docelową ${p}$ , jest określone poniższym wzorem:

$\qquad\begin{split}\chi^{2}=\chi^{{2}}({q},{p})&=\int\frac{(q(x)-{p}(x))^{{2}}}{q(x)}{\rm d}x\\ &=\mathbb{E}_{{q}}\left(\frac{q(X)-{p}(X)}{q(X)}\right)^{{2}}=\mathbb{E}_{{q}}\left(1-w(X)\right)^{{2}}\\ &={\rm Var}_{{q}}w(X).\end{split}$

Niestety, to ,,odchylenie” nie jest odległością, bo nie spełnia warunku symetrii. Niemniej widać, że małe wartości ${\rm Var}_{{q}}w(X)$ świadczą o bliskości ${q}$ i ${p}$ . Zauważmy, jeszcze, że ${\rm Var}_{{q}}w(X)$ jest wariancją asymptotyczną estymatora IS1 dla funkcji $f(x)=1$ .

Zgodnym estymatorem wielkości $\chi^{2}$ jest

$\qquad\hat{\chi}^{2}=\frac{1}{n-1}\sum _{{i=1}}^{{n}}\frac{\left(\bar{W}_{{n}}-W_{{i}}\right)^{{2}}}{\phantom{\big(}\bar{W}^{{2}}\phantom{\big(}}.$

Dzielenie przez $n-1$ (a nie $n$ ) wynika z tradycji. Dzielenie przez $\bar{W}^{{2}}$ powoduje, że estymator nie wymaga znajomości czynnika normującego gęstość ${p}$ . Jest to w istocie estymator typu IS2. Estymator typu IS1 jest równy po prostu $\sum(1-W_{{i}})^{2}/n$ , ale można go stosować tylko gdy znamy czynnik normujący.

Wprowadźmy teraz bardziej formalnie pojęcie ważonej próbki. Niech $(X,W)$ będzie parą zmiennych losowych o wartościach w ${\cal X}\times[0,\infty[$ . Jeśli rozkład brzegowy $X$ ma gęstość $q$ zaś docelowa gęstość jest postaci $p(x)=p^{\prime}(x)/Z$ to żądamy aby

$\qquad\mathbb{E}(W|X=x)=\frac{p^{\prime}(x)}{q(x)}$

dla prawie wszystkich $x$ . Mówimy wtedy, że próbka jest dopasowana do rozkładu $p$ . Idea jest taka jak w losowaniu istotnym IS2: dopasowanie gwarantuje, że dla każdej funkcji $f$ mamy $\mathbb{E}_{q}f(X)W=Z\mathbb{E}_{p}f(X)$ . Nie jest konieczne, aby $W$ było deterministyczną funkcją $X$ . Pozostawiamy sobie możliwość generowania wag losowo.

Stała normująca $Z$ jest na ogół nieznana, a więc waga $W$ jest określona z dokładnością do proporcjonalności. W poprzednich podrozdziałach dla podkreślenia opatrywaliśmy takie nieunormowane wagi znaczkiem ,,prim” ale teraz to pomijamy. Gdy generujemy ciąg ważonych zmiennych losowych $(X_{{i}},W_{{i}}),\ldots,(X_{{n}},W_{{n}})$ , to oczywiście wymaga się aby $\mathbb{E}(W_{{i}}|X_{{i}}=x)=Zp(x)/q(x)$ , gdzie stała $Z$ musi być jednakowa dla wszystkich $i=1,\ldots,n$ . Jeśli założy się niezależność par $(X_{{i}},W_{{i}})$ , to zgodność i asymptotyczną normalność estymatora $\hat{\theta}_{n}^{{\rm IS2}}$ danego równaniem (8.2) wykazuje się tak samo jak poprzednio. Asymptotyczna wariancja jest równa $Z^{{-2}}{\rm Var}_{{q}}W(f(X)-\theta)$ ; czynnik $Z^{{-2}}$ pojawia się dlatego, że operujemy nieunormowanymi wagami.

Jest jednak dużo ciekawych zastosowań, w których punkty $X_{{i}}$ są zależne. Wtedy nawet zgodność estymatora $\hat{\theta}_{n}^{{\rm IS2}}$ nie jest automatyczna. Asymptotyczna normalność może zachodzić z graniczną wariancją zupełnie inną niż w przypadku niezależnym lub nie zachodzić wcale.

8.3. Ważona eliminacja

Interesujące jest powiązanie idei eliminacji z ważeniem próbek. Czysta metoda eliminacji pracuje w sytuacji, gdy gęstość ,,instrumentalna” majoryzuje funkcję proporcjonalną do gęstości docelowej, czyli $Zp(x)=p^{\prime}(x)\leqslant q(x)$ . Jeśli ten warunek nie jest spełniony, to można naprawić odpowiednio ważąc wylosowaną próbkę. Dokładniej, algorytm ważonej eliminacji (Rejection Control, w skrócie RC) jest następujący.

repeat

Gen $Y\sim q$ ;

$W:=p^{\prime}(Y)/q(Y)$ ;

if $W\leqslant 1$ then

begin

$W:=1$ ;

Gen $U\sim{\rm U}(0,1)$

end;

until ( $W>1$ or $U<W$ )

$X:=Y$

Dowód poprawności algorytmu

Zauważmy, że prawdopodobieństwo akceptacji $Y$ w tym algorytmie jest równe $a(Y)$ , gdzie

$\qquad a(x)=\frac{p^{\prime}(x)\land q(x)}{q(x)}=\begin{cases}\dfrac{p^{\prime}(x)}{q(x)}&\textrm{ jeśli }p^{\prime}(x)\leqslant q(x);\\ 1&\textrm{ jeśli }p^{\prime}(x)>q(x).\end{cases}$

Rozkład $X$ na wyjściu algorytmu ma zatem gęstość $\propto p^{\prime}(x)\land q(x)$ . Waga $W=w(Y)$ jest więc ,,dopasowana” do rozkładu $p$ w sensie zdefiniowanym w poprzednim podrozdziale, bo

$\qquad w(x)=\frac{p^{\prime}(x)}{p^{\prime}(x)\land q(x)}=\begin{cases}1&\textrm{ jeśli }p^{\prime}(x)\leqslant q(x);\\ \dfrac{p^{\prime}(x)}{q(x)}&\textrm{ jeśli }p^{\prime}(x)>q(x).\end{cases}$

∎

Zauważmy jeszcze, jak należy modyfikować wagi w RC, jeśli na wejściu mamy próbkę ważoną, dopasowaną do $p$ . Jeżeli punkt $Y$ , pochodzący z rozkładu $q$ , ma na wejściu wagę $W_{{Y}}=p^{\prime}(Y)/q(Y)$ , to na wyjściu zaakceptowany punkt $X=Y$ ma rozkład $\propto p^{\prime}(x)\land q(x)$ i powinien mieć wagę $W_{{X}}=p^{\prime}(Y)/(p^{\prime}(Y)\land q(Y))$ . Widać stąd, że

$\qquad W_{Y}=\frac{W_{X}}{a(Y)},$

gdzie $a(\cdot)$ jest napisanym wyżej prawdopodobieństwem akceptacji w RC.

Przykład 8.1 (Nie-samo-przecinające się błądzenia)

Po angielsku nazywają się Self Avoiding Walks, w skrócie SAW. Niech $\mathbb{Z}^{{d}}$ będzie $d$ -wymiarową kratą całkowitoliczbową. Mówimy, że ciąg $s=(0=s_{{0}},s_{{1}},\ldots,s_{{k}})$ punktów kraty jest SAW-em jeśli

każde dwa kolejne punkty $s_{{i-1}}$ i $s_{{i}}$ sąsiadują ze sobą, czyli różnią się o $\pm 1$ na dokładnie jednej współrzędnej,
żadne dwa punkty nie zajmują tego samego miejsca, czyli $s_{{i}}\not=s_{{j}}$ dla $i\not=j$ .

Zbiór wszystkich SAW-ów o $k$ ogniwach w $\mathbb{Z}^{{d}}$ oznaczymy $\text{SAW}_{{k}}^{{d}}$ , a dla $d$ i $k$ ustalonych w skrócie SAW. Przykład $s\in\text{SAW}_{{15}}^{{2}}$ widać na rysunku.

Rys. 8.1. Przykład SAW-a.

Natychmiast nasuwa się bardzo proste pytanie:

Jak policzyć SAW-y, czyli obliczyć $L=L_{{k,d}}=|\text{SAW}_{{k}}^{{d}}|$ ?

Zainteresujmy się teraz ,,losowo wybranym SAW-em”. Rozumiemy przez to zmienną losową $S$ o rozkładzie jednostajnym w zbiorze $\text{SAW}_{{k}}^{{d}}$ , czyli taką, że $\mathbb{P}(S=s)=1/L_{{k,d}}$ dla każdego $s\in\text{SAW}_{{k}}^{{d}}$ . W skrócie, $S\sim{\rm U}(\textrm{SAW})$ . Niech $\delta(s_{{k}})$ oznacza odległość euklidesową końca SAW-a, czyli punktu $s_{{k}}$ , od początku, czyli punktu $0$ . Na przykład dla lańcuszka widocznego na rysunku mamy $\delta(s_{{15}})=\sqrt{5}$ . Można zadać sobie pytanie, jaka jest średni kwadrat takiej odległości, czyli

Jak obliczyć $\overline{\delta^{2}}=\overline{\delta _{{d,k}}^{2}}=\mathbb{E}\delta(S_{{k}})^{{2}}$ , przy założeniu, że $S\sim{\rm U}(\textrm{SAW})$ ?

Można zastosować prostą metodę MC, czyli eliminację. Niech $\text{WALK}_{{k}}^{{d}}$ oznacza zbiór wszystkich ,,błądzeń”, czyli ciągów $s=(0=s_{{0}},s_{{1}},\ldots,s_{{k}})$ niekonieczne spełniających warunek ,, $s_{{i}}\not=s_{{j}}$ dla $i\not=j$ ”. Oczywiście $|\text{WALK}_{{k}}^{{d}}|=(2d)^{{k}}$ i metoda generowania ,,losowego błądzenia” (z rozkładu ${\rm U}(\text{WALK})$ ) jest bardzo prosta: kolejno losujemy pojedyncze kroki, wybierając zawsze jedną z $2d$ możliwości. Żeby otrzymać ,,losowy SAW”, stosujemy eliminację. Ten sposób pozwala w zasadzie estymować $\overline{\delta^{2}}$ (przez uśrednianie długości zaakceptowanych błądzeń) oraz $\text{SAW}/\text{WALK}$ (przez zanotowanie frakcji akceptowanych błądzeń). Niestety, metoda jest bardzo nieefektywna, bo dla dużych $k$ prawdopodobieństwo akceptacji szybko zbliża sie do zera.

Rys. 8.2. Generowanie SAW-a metodą ,,wzrostu”.

Metoda ,,wzrostu” zaproponowana przez Rosenbluthów polega na losowaniu kolejnych kroków błądzenia spośród ,,dopuszczalnych punktów”, to znaczy punktów wcześniej nie odwiedzonych. W każdym kroku, z wyjątkiem pierwszego mamy co najwyżej $2d-1$ możliwości. W błądzeniu widocznym na rysunku kolejne kroki wybieraliśmy spośród:

$4,3,3,3,\quad 2,3,2,2,\quad 3,2,3,3,\quad 2,1,3,2$

możliwych. Nasz SAW został zatem wylosowany z prawdopodobieństwem

$\frac{1}{4}\cdot\frac{1}{3}\cdot\frac{1}{3}\cdot\frac{1}{3}\cdot\quad\frac{1}{2}\cdot\frac{1}{3}\cdot\frac{1}{2}\cdot\frac{1}{2}\cdot\quad\frac{1}{3}\cdot\frac{1}{2}\cdot\frac{1}{3}\cdot\frac{1}{3}\cdot\quad\frac{1}{2}\cdot\frac{1}{1}\cdot\frac{1}{3}\cdot\frac{1}{2}$

Powiedzmy ogólniej, że przy budowaniu SAW-a $s=(0=s_{{0}},s_{{1}},\ldots,s_{{k}})$ mamy kolejno

$n_{1}=2d,n_{2},\ldots,n_{k}$

możliwości (nie jest przy tym wykluczone, że w pewnym kroku nie mamy żadnej możliwości, $n_{i}=0$ ). Używając terminologii i oznaczeń związanych z losowaniem istotnym powiemy, że

$q(s)=\frac{1}{n_{1}}\cdot\frac{1}{n_{2}},\cdots\frac{1}{n_{k}}.$

jest gęstością instrumentalną dla $s\in\text{SAW}$ , gęstość docelowa jest stała, równa $p(s)=1/Z$ , zatem funkcja podziału jest po prostu liczbą SAW-ów: $Z=|\text{SAW}|$ . Wagi przypisujemy zgodnie ze wzorem $w(s)={n_{1}}\cdot{n_{2}}\cdots{n_{k}}$ (jeśli wygenerowanie SAW-a się nie udało, $n_{i}=0$ dla pewnego $i$ , to waga jest zero). Niech teraz $S(1),\ldots,S(n)$ będą niezależnymi błądzeniami losowanymi metodą Rosenbluthów. Zgodnie z ogólnymi zasadami losowania istotnego,

$\frac{\sum w(S(i))\delta(S(i))}{\sum w(S(i))}\quad\text{jest estymatorem }\overline{\delta^{2}},$

$\sum w(S(i))\quad\text{jest estymatorem }|SAW|.$

Należy przy tym uwzględniać w tych wzorach błądzenia o wadze zero, czyli ,,nieudane SAW-y”.

Przykład 8.2 (Prawdopodobieństwo ruiny i wykładnicza zamiana miary)

Rozpatrzymy najprostszy model procesu opisującego straty i przychody w ubezpieczeniowej ,,teorii ryzyka”. Niech $Y_{{1}},Y_{{2}},\ldots$ będą zmiennymi losowymi oznaczającymi straty netto (straty $-$ przychody) towarzystwa ubezpieczeniowego w kolejnych okresach czasu. Założymy (co jest dużym uproszczeniem), że te zmienne są niezależne i mają jednakowy rozkład o gęstości $p(y)$ . Tak zwana ,,nadwyżka ubezpieczyciela” na koniec $n$ -go roku jest równa

$u-S_{n}=u-\sum _{{i=1}}^{n}Y_{i},$

gdzie $u$ jest rezerwą początkową. Interesuje nas prawdopodobieństwo zdarzenia, polegającego na tym, że $u-S_{n}<0$ dla pewnego $n$ . Mówimy wtedy (znowu w dużym uproszczeniu) o ,,ruinie ubezpieczyciela”. Wygodnie jest przyjąć następujące oznaczenia i konwencje. Zmienna losowa

$R=\begin{cases}\min\{ n:S_{n}>u\}&\text{ jeśli takie $n$ istnieje};\\ \infty&\text{ jeśli $S_{n}\leqslant u$ dla wszystkich $n$}\end{cases}$

oznacza czas oczekiwania na ruinę, przy czym jeśli ruina nigdy nie nastąpi to ten czas uznajemy za nieskończony. Przy takiej umowie, prawdopodobieństwo ruiny możemy zapisać jako $\psi=\mathbb{P}_{p}(R<\infty)$ . Wskaźnik $p$ przy symbolu wartości oczekiwanej przypomina, że chodzi tu o ,,oryginalny” proces, dla którego $Y_{{i}}\sim p$ .

Obliczenie $\psi$ analitycznie jest możliwe tylko w bardzo specjalnych przykładach. Motody numeryczne istnieją, ale też nie są łatwe. Pokażemy sposób obliczania $\psi$ metodą Monte Carlo, który stanowi jeden z najpiękniejszych, klasycznych, przykładów losowania istotnego. Przyjmiemy bardzo rozsądne założenie, że $\mathbb{E}_{p}Y_{{i}}<0$ . Funkcję tworzącą momenty, która odpowiada gęstości $p$ określamy wzorem

${M}_{p}(t)=\mathbb{E}_{p}{\rm e}^{{tY_{{i}}}}=\int _{{-\infty}}^{\infty}{\rm e}^{{ty}}p(y){\rm d}y.$

Założymy, że ta funkcja przyjmuje wartości skończone przynajmniej w pewnym otoczeniu zera i istnieje takie $r>0$ , że

${M}_{p}(t)=1.$

Liczba $r$ jest nazywana współczynnikiem dopasowania i odgrywa tu kluczową rolę.

Metoda wykładniczej zamiany miary jest specjalnym przypadkiem losowania istotnego. Zeby określić instrumentalny rozkład prawdopodobieństwa, połóżmy

$q(y)={\rm e}^{{ry}}p(y).$

Z definicji współczynnika dopasowania wynika, że $q$ jest gęstością prawdopodobieństwa, to znaczy $\int q(y){\rm d}y=1$ . Generuje się ciąg $Y_{1},Y_{2},\ldots$ jednakowo rozłożonych zmiennych losowych o gęstości $q$ . Dla utworzonego w ten sposób ,,instrumentalnego procesu” używać będziemy symboli $\mathbb{P}_{q}$ i $\mathbb{E}_{q}$ . Zauważmy, że

$\begin{split}\mathbb{E}_{q}Y_{i}=\int yq(y){\rm d}y&=\int y{\rm e}^{{ry}}p(y){\rm d}y\\ &=\frac{{\rm d}}{{\rm d}r}\int{\rm e}^{{ry}}p(y){\rm d}y\\ &={M}_{p}^{\prime}(r)>0,\end{split}$

ponieważ funkcja tworząca momenty ${M}_{p}(\cdot)$ jest wypukła, ${M}_{p}^{\prime}(0)<0$ i ${M}_{p}(0)={M}_{p}(r)=1$ . Po zamianie miary, proces $u-S_{n}$ na mocy Prawa Wielkich Liczb zmierza prawie na pewno do minus nieskończoności, a zatem ruina następuje z prawdopodobieństwem jeden, $\mathbb{P}_{q}(R<\infty)=1$ . Pokażemy, jak wyrazić prstwo ruiny dla oryginalnego procesu w terminach procesu instrumentalnego. Niech

${\cal R}_{n}=\{(y_{1},\ldots,y_{n}):y_{1}\leqslant u,\ldots,y_{1}+\cdots+y_{{n-1}}\leqslant u,y_{1}+\cdots+y_{{n-1}}+y_{n}>u\},$

Innymi słowy, zdarzenie $\{ R=n\}$ zachodzi gdy $(Y_{1},\ldots,Y_{n})\in{\cal R}_{n}$ . Mamy zatem

$\begin{split}\mathbb{P}_{p}(R=n)&={\int\cdots\int}_{{{\cal R}_{n}}}p(y_{1})\cdots p(y_{n}){\rm d}y_{1}\cdots{\rm d}y_{n}\\ &={\int\cdots\int}_{{{\cal R}_{n}}}{\rm e}^{{-ry_{1}}}q(y_{1})\cdots{\rm e}^{{-ry_{n}}}q(y_{n}){\rm d}y_{1}\cdots{\rm d}y_{n}\\ &={\int\cdots\int}_{{{\cal R}_{n}}}{\rm e}^{{-r(y_{1}+\cdots+y_{n})}}q(y_{1})\cdots q(y_{n}){\rm d}y_{1}\cdots{\rm d}y_{n}\\ &=\mathbb{E}_{q}{\rm e}^{{rS_{n}}}\mathbb{I}(R=n).\end{split}$

Weźmy sumę powyższych równości dla $n=1,2,\ldots$ i skorzystajmy z faktu, że $\sum _{{n=1}}^{\infty}\mathbb{P}_{q}(R=n)=1$ . Dochodzimy do wzoru

$\qquad\mathbb{P}_{p}(R<\infty)=\mathbb{E}_{q}\exp\{-rS_{{R}}\}={\rm e}^{{-ru}}\mathbb{E}_{q}\exp\{-r(S_{{R}}-u)\}.$

(8.3)

Ten fakt jest podstawą algorytmu Monte Carlo:

$\hat{\psi}:=0$ ; [ $\hat{\psi}$ będzie estymatorem prawdopodobieństwa ruiny ]

for $m:=1$ to $k$ do

begin

$n:=0$ ; $S:=0$ ;

repeat

Gen $Y\sim q$ ; $S:=S+Y$ ;

until $S>u$ ;

$Z:=\exp\{-r(S-u)\}$ ; $\hat{\psi}:=\hat{\psi}+Z$

end

$\hat{\psi}:={\rm e}^{{-ru}}\hat{\psi}/k$

Algorytm jest prosty i efektywny. Trochę to zadziwiające, że w celu obliczenia prawdopodobieństwa ruiny generuje się proces dla którego ruina jest pewna. Po chwili zastanowienia można jednak zauważyć, że wykładnicza zamiana miary realizuje podstawową ideę losowania istotnego: rozkład instrumentalny ,,naśladuje” proces docelowy ograniczony do zdarzenia $\{ R<\infty\}$ .

Ciekawe, że wykładnicza zamiana miary nie tylko jest techniką Monte Carlo, ale jest też techniką dowodzenia twierdzeń! Aby się o tym przekonać, zauważmy, że ,,po drodze” udowodniliśmy nierówność $\psi\leqslant{\rm e}^{{-ru}}$ (wynika to z podstawowego wzoru (8.3) gdyż $\exp\{-r(S_{R}-u)\}\leqslant 1$ ). Jest to sławna nierówność Lundberga i wcale nie jest ona oczywista.