Zagadnienia

14.1 Wstęp, metody niedeterministyczne
14.2 Klasyczna metoda Monte Carlo
- 14.2.1 Definicja i błąd
- 14.2.2 Całkowanie z wagą
14.3 Redukcja wariancji
- 14.3.1 Losowanie warstwowe
- 14.3.2 Funkcje kontrolne
14.4 Generowanie liczb (pseudo-)losowych

14. Metody Monte Carlo

14.1. Wstęp, metody niedeterministyczne

Poprzedni wykład zakończyliśmy pesymistycznym twierdzeniem 13.4, że nie istnieją efektywne metody numerycznego całkowania funkcji wielu zmiennych, ponieważ ma miejsce zjawisko przekleństwa wymiaru. Zwróćmy jednak uwagę na to, że fakt istnienia przekleństwa wymiaru stwierdziliśmy przy założeniach, że:

(i) model obliczeniowy jest deterministyczny,
(ii) funkcje podcałkowe są $r$ -krotnie różniczkowalne po każdej zmiennej.

Można mieć nadzieję, że przekleństwo wymiaru zniknie, albo zostanie złagodzone, gdy przynajmniej jedno z tych założeń nie będzie spełnione.

Ten wykład poświęcimy (klasycznej) metodzie Monte Carlo numerycznego całkowania, która jest przykładem metody niedeterministycznej, tzn. takiej, która oblicza wynik wykorzystując zjawiska losowe. Chociaż może to brzmieć dziwnie, to właśnie niedeterministyczne zachowanie metody pozwala pokonać przekleństwo wymiaru.

Opisana dalej klasyczna metoda Monte Carlo związana jest ściśle ze Stanisławem Ulamem, uczniem Stefana Banacha i reprezentantem Lwowskiej Szkoły Matematycznej. Ulam zastosował metodę Monte Carlo do obliczania skomplikowanych całek w ramach ,,Manhattan Project” w Los Alamos (USA), w czasie II Wojny światowej.

14.2. Klasyczna metoda Monte Carlo

14.2.1. Definicja i błąd

Tak jak w poprzednim rozdziale chcemy obliczyć całkę

$I_{d}(f)\,:=\,\int _{{[0,1]^{d}}}f(\vec{x})\, d\vec{x}\,=\,\underbrace{\int _{0}^{1}\int _{0}^{1}\cdots\int _{0}^{1}}_{d}f(x_{1},x_{2},\ldots,x_{d})\, dx_{1}dx_{2}\cdots dx_{d}.$

Zakładamy przy tym, że $f:[0,1]^{d}\to\mathbb{R}$ jest funkcją, której kwadrat jest całkowalny,

$\int _{{[0,1]^{d}}}|f(\vec{x})|^{2}\, d\vec{x}\,<\,\infty.$

Definicja 14.1

Klasyczna metoda Monte Carlo polega na przybliżeniu $I_{d}(f)$ średnią arytmetyczną wartości funkcji $f$ w losowo wybranych punktach, tzn.

$MC_{{d,N}}(f)\,=\, MC_{{d,N}}(f;\vec{t}_{1},\vec{t}_{2},\ldots,\vec{t}_{N})\,:=\,\frac{1}{N}\cdot\sum _{{j=1}}^{N}f(\vec{t}_{j}),$

gdzie $\vec{t}_{1},\vec{t}_{2},\ldots,\vec{t}_{N}$ są punktami wylosowanymi niezależnie od siebie, każdy zgodnie z rozkładem jednostajnym na $[0,1]^{d}$ .

Konsekwencją zastosowania losowości jest to, że przy różnych realizacjach metody otrzymujemy różne wyniki, w zależności od wyboru punktów $\vec{t}_{j}$ . Wynik $MC_{{d,N}}(f)$ jest więc zmienną losową, której wartość oczekiwana wynosi

	$\displaystyle E\left(MC_{{d,N}}(f)\right)$	$\displaystyle=$	$\displaystyle\int _{{[0,1]^{{d\cdot N}}}}MC_{{d,N}}(f;\vec{t}_{1},\ldots,\vec{t}_{N})\, d\vec{t}_{1}\cdots d\vec{t}_{N}$
		$\displaystyle=$	$\displaystyle\frac{1}{N}\sum _{{j=1}}^{N}\int _{{[0,1]^{d}}}f(\vec{t})\, d\vec{t}\,=\, I_{d}(f).$

Ponieważ różnica $I_{d}(f)-MC_{{d,N}}(f)$ jest też zmienną losową, za błąd metody Monte Carlo dla danej funkcji $f$ przyjmiemy odchylenie standardowe,

$e(f;MC_{{d,N}})\,:=\,\sqrt{E\left(I_{d}(f)-MC_{{d,N}}(f)\right)^{2}}.$

Twierdzenie 14.1

Dla danej funkcji $f$ błąd metody Monte Carlo wynosi

$e(f;MC_{{d,N}})\,=\,\frac{\sigma(f)}{\sqrt{N}},$

gdzie

$\sigma(f)\,:=\,\sqrt{I_{d}(f^{2})-I_{d}^{2}(f)}$

jest wariancją funkcji $f$ .

Zanim przystąpimy do dowodu zauważmy, że $\sigma(f)$ jest dobrze zdefiniowaną wielkością, bowiem nierówność

$|I_{d}(f)|\le\sqrt{I_{d}(f^{2})}$

jest szczególnym przypadkiem znanej nierówności Schwarza dla całek.

Dowód

Oznaczmy, dla uproszczenia, zmienną losową $X=MC_{{d,N}}(f)$ . Wtedy

$E(X-E(X))^{2}\,=\, E\left(X(X-E(X))-E(X)(X-E(X))\right)\,=\, E(X^{2})-E^{2}(X).$

(14.1)

Ponadto

	$\displaystyle E(X^{2})$	$\displaystyle=$	$\displaystyle E\left(\frac{1}{N}\left(\sum _{{j=1}}^{N}f(\vec{t}_{j})\right)^{2}\right)\,=\,\frac{1}{N^{2}}E\left(\sum _{{j=1}}^{N}f^{2}(\vec{t}_{j})+\sum _{{i\ne j}}f(\vec{t}_{i})f(\vec{t}_{j})\right)$
		$\displaystyle=$	$\displaystyle\frac{1}{N^{2}}\left(NI_{d}(f^{2})+(N^{2}-N)I_{d}^{2}(f)\right)\,=\,\frac{1}{N}I_{d}(f^{2})+\left(1-\frac{1}{N}\right)I_{d}^{2}(f),$

gdzie skorzystaliśmy z niezależności zmiennych losowych $f(t_{i})$ i $f(t_{j})$ dla $i\ne j$ . Stąd i z (14.1) dostajemy

$e^{2}(f;MC_{{d,N}})\,=\, E(X-E(X))^{2}\,=\,\frac{1}{N}I_{d}(f^{2})+\left(1-\frac{1}{N}\right)I_{d}^{2}(f)-I_{d}^{2}(f)\,=\,\frac{1}{N}\left(I_{d}(f^{2})-I_{d}^{2}(f)\right),$

co kończy dowód.

∎

Uwaga 14.1

Zauważmy, że w dowodzie pokazaliśmy przy okazji nierówność Schwarza posługując się narzędziami rachunku prawdopodobieństwa.

Twierdzenie (14.1) mówi, że błąd metody Monte Carlo jest proporcjonalny do $N^{{-1/2}}$ przy bardzo słabych wstępnych założeniach na funkcję (jedynie całkowalność kwadratu funkcji). Jest to istotna poprawa w porównaniu do błędu $N^{{-r/d}}$ dla metod deterministycznych. W szczególności ważne jest, że wykładnik $1/2$ przy $N^{{-1}}$ jest niezależny od wymiaru $d$ , a konsekwencją tego pokonanie przekleństwa wymiaru.

Dziwnym może wydawać się, że przekleństwo wymiaru można zlikwidować używając metod niedeterministycznych (losowych). Jednak niczego nie ma za darmo. Należy pamiętać, że jest to możliwe za cenę niepewności wyniku. O ile bowiem metoda deterministyczna produkuje zawsze ten sam wynik, metoda niedeterministyczna (taka jak Monte Carlo) produkuje różne wyniki zależnie od konkretnych realizacji zmiennych losowych. Dlatego, mimo iż błąd oczekiwany jest proporcjonalny do $N^{{-1/2}}$ to nie mamy całkowitej pewności, że przy konkretnej realizacji otrzymany wynik jest tego samego rzędu. Z tego punktu widzenia warto przytoczyć następującą równość, która wynika z centralnego twierdzenia granicznego; mianowicie, dla dowolnych $c_{1}<c_{2}$ mamy

$\lim _{{N\to\infty}}\mbox{Prob}\left(\frac{c_{1}\sigma(f)}{\sqrt{N}}\le I_{d}(f)-MC_{{d,N}}(f)\le\frac{c_{2}\sigma(f)}{\sqrt{N}}\right)\,=\,\frac{1}{\sqrt{2\pi}}\int _{{c_{1}}}^{{c_{2}}}e^{{-t^{2}/2}}dt,$

gdzie Prob oznacza prawdopodobieństwo względem rozkładu jednostajnego na $[0,1]^{{dN}}$ .

14.2.2. Całkowanie z wagą

Deterministyczne metody interpolacyjne z poprzedniego rozdziału można stosować jedynie do całkowania na $d$ -wymiarowych prostokątach. Metoda Monte Carlo ma oprócz wymienionych również i tą zaletę, że łatwo ją uogólnić na przypadek całkowania z wagą. Dla przybliżenia wartości

$I_{{d,\omega}}(f)\,:=\,\int _{{\mathbb{R}^{d}}}f(\vec{x})\omega(\vec{x})\, d\vec{x},\qquad\mbox{gdzie}\qquad\int _{{\mathbb{R}^{d}}}\omega(\vec{x})\, d\vec{x}\,=:\, W\,<\,\infty,$

możemy bowiem zastosować wzór

$MC_{{d,N}}^{\omega}(f)\,:=\,\frac{W}{N}\cdot\sum _{{j=1}}^{N}f(\vec{t}_{j}),$

przy czym $\vec{t}_{1},\ldots,\vec{t}_{N}$ są tym razem punktami wybranymi losowo i niezależnie od siebie, zgodnie z rozkładem na $\mathbb{R}^{d}$ o gęstości $\omega/W$ .

Adaptując odpowiednio dowód twierdzenia 14.1 otrzymujemy następujące wyrażenie na błąd uogólnionej metody Monte Carlo.

Twierdzenie 14.2

Niech $I_{{d,\omega}}(f^{2})=\int _{{\mathbb{R}^{d}}}f^{2}(\vec{x})\omega(\vec{x})\, d\vec{x}<\infty$ . Wtedy

$e(f;MC_{{d,N}}^{\omega})\,=\,\frac{\sigma _{\omega}(f)}{\sqrt{N}},$

gdzie

$\sigma _{\omega}(f)\,=\,\sqrt{\, W\, I_{{d,\omega}}(f^{2})-I^{2}_{{d,\omega}}(f)}.$

14.3. Redukcja wariancji

Zauważyliśmy, że zaletą metody Monte Carlo jest nie tylko jej prostota, ale również to, że błąd średni wynosi $\sigma _{\omega}(f)N^{{-1/2}}$ . Naturalnym jest teraz pytanie, czy błędu tego nie można poprawić. Temu celowi służą metody redukcji wariancji, które polegają w ogólności na redukcji czynnika $\sigma _{\omega}(f)$ . Spośród wielu technik redukcji wariancji skupimy uwagę na dwóch: losowaniu warstwowemu oraz funkcjach kontrolnych. Dla uproszczenia będziemy zakładać, że całkujemy z wagą jednostkową na kostce

$D=[0,1]^{d}.$

14.3.1. Losowanie warstwowe

Podzielmy obszar całkowania $D$ na $K$ rozłącznych podzbiorów $D_{i}$ tak, że

$D=\bigcup _{{i=1}}^{K}D_{i}$

i zastosujmy Monte Carlo do całkowania po każdym $D_{i}$ , tzn. całkę $\int _{D}f(\vec{x})\, d\vec{x}$ przybliżymy wielkością

$\overline{M}C_{{d,N}}(f)\,:=\,\sum _{{i=1}}^{K}MC_{{d,N_{i}}}^{{(i)}}(f),$

gdzie $MC_{{d,N_{i}}}^{{(i)}}$ jest metodą Monte Carlo zastosowaną do całki

$I_{d}^{{(i)}}(f)\,:=\,\int _{{D_{i}}}f(\vec{x})\, d\vec{x},\qquad 1\le i\le K,$

oraz $N=\sum _{{i=1}}^{K}N_{i}$ .

Oznaczmy przez $|D_{i}|$ objętość $d$ -wymiarową podzbioru $D_{i}$ . Ponieważ zmienne losowe $I_{d}^{{(i)}}(f)-MC_{{d,N_{i}}}^{{(i)}}(f)$ są parami niezależne dla $1\le i\le K$ , na podstawie twierdzenia 14.2 mamy

$\displaystyle E\left((I_{d}(f)-\overline{M}C_{{d,N}}(f))^{2}\right)$	$\displaystyle=$	$\displaystyle E\left(\left(\sum _{{i=1}}^{K}I_{d}^{{(i)}}(f)-MC_{{d,N_{i}}}(f)\right)^{2}\right)$
	$\displaystyle=$	$\displaystyle\sum _{{i=1}}^{K}E\left((I_{d}^{{(i)}}(f)-MC_{{d,N_{i}}}^{{(i)}}(f))^{2}\right)$
	$\displaystyle=$	$\displaystyle\sum _{{i=1}}^{K}\frac{1}{N_{i}}\left(\|D_{i}\|\, I_{d}^{{(i)}}(f^{2})-(I_{d}^{{(i)}}(f))^{2}\right).$

Przyjmijmy teraz, że

$N_{i}\,=\,|D_{i}|\cdot N,\qquad 1\le i\le K,$

przy czym dla uproszczenia (ale bez utraty ogólności) zakładamy, że wielkości te są całkowite. Wtedy otrzymujemy

$e(f,\overline{M}C_{{d,N}})\,=\,\frac{1}{\sqrt{N}}\cdot\sqrt{I_{d}(f^{2})-\sum _{{i=1}}^{K}\frac{1}{|D_{i}|}(I_{d}^{{(i)}}(f))^{2}}.$

(14.2)

Błąd tak zdefiniowanej metody $\overline{M}C_{{d,N}}$ nie jest większy od błędu klasycznej metody $MC_{{d,N}}$ z Twierdzenia 14.1.

Twierdzenie 14.3

Dla dowolnej funkcji $f$ takiej, że $I_{d}(f^{2})<\infty$ mamy

$e(f,\overline{M}C_{{d,N}})\,\le\, e(f,MC_{{d,N}}),$

przy czym równość zachodzi tylko wtedy gdy iloraz $I_{d}^{{(i)}}(f)/|D_{i}|$ jest stały, niezależnie od $i$ .

Dowód

Rzeczywiście, oznaczając

$a_{i}=\sqrt{|D_{i}|},\qquad b_{i}=\frac{I_{d}^{{(i)}}(f)}{\sqrt{|D_{i}|}},$

oraz wykorzystując nierówność Schwarza dla ciągów mamy

$I_{d}^{2}(f)\,=\,\left(\sum _{{i=1}}^{K}a_{i}b_{i}\right)^{2}\,\le\,\left(\sum _{{i=1}}^{K}a_{i}^{2}\right)\left(\sum _{{i=1}}^{K}b_{i}^{2}\right)\,=\,\sum _{{i=1}}^{K}b_{i}^{2}\,=\,\sum _{{i=1}}^{K}\frac{1}{|D_{i}|}(I_{d}^{{(i)}}(f))^{2},$

przy czym równość zachodzi tylko wtedy gdy wektory $(a_{1},\ldots,a_{K})^{T}$ i $(b_{1},\ldots,b_{K})^{T}$ są liniowo zależne, co jest równoważne warunkowi w treści twierdzenia.

Prawdziwość tezy pokazuje teraz porównanie wzorów na błędy obu metod.

∎

Widzimy, że stosując losowanie warstwowe z ustalonym podziałem na $K$ podzbiory $D_{i}$ możemy co prawda zmiejszyć błąd, ale szybkość zbieżności $N^{{-1/2}}$ pozostaje ta sama. A czy można poprawić zbieżność stosując różne podziały dla różnych wartości $N$ ? Okazuje się, że tak, o ile założymy pewną regularność funkcji $f$ .

Aby to uzyskać, najpierw przekształcimy wzór (14.2) na błąd metody $\overline{M}C_{{d,N}}$ do postaci

$e(f,\overline{M}C_{{d,N}})\,=\,\frac{1}{\sqrt{N}}\sqrt{\sum _{{i=1}}^{K}I_{d}^{{(i)}}((f-c_{i})^{2})}.$

(14.3)

gdzie

$c_{i}:=\frac{I_{d}^{{(i)}}(f)}{|D_{i}|},\qquad 1\le i\le K.$

Załóżmy teraz, że $f$ spełnia warunek Lipschitza ze stałą $L$ ,

$|f(\vec{x})-f(\vec{y})|\,\le\, L\cdot\|\vec{x}-\vec{y}\| _{\infty},\qquad\vec{x},\vec{y}\in D.$

Wtedy istnieją $\vec{x}_{i}\in\overline{D}_{i}$ takie, że $f(\vec{x}_{i})=c_{i}$ , a stąd i z lipschitzowskości $f$ mamy, że dla dowolnego $\vec{x}\in D_{i}$

$|f(\vec{x})-c_{i}|\,\le\, L\cdot\|\vec{x}-\vec{x}_{i}\| _{\infty}\,\le\, L\cdot{\rm diam}_{\infty}(D_{i}),$

gdzie

${\rm diam}_{\infty}(D_{i})\,:=\,\sup\left\{\|\vec{x}-\vec{y}\| _{\infty}\,:\;\vec{x},\vec{y}\in D_{i}\right\}$

jest średnicą zbioru $D_{i}$ w normie max. W konsekwencji, ze wzoru (14.3) dostajemy następujące oszacowanie błędu:

$e(f,\overline{M}C_{{d,N}})\,\le\,\frac{L}{\sqrt{N}}\sqrt{\sum _{{i=1}}^{K}|D_{i}|\,{\rm diam}^{2}(D_{i})}.$

Ustalmy teraz równomierny podział kostki $D$ na $K=N$ podkostek $D_{i}$ , każda o krawędzi długości $N^{{-1/d}}$ (zakładamy, bez zmniejszenia ogólności, że $N^{{1/d}}$ jest całkowita) tak, że nasza metoda aproksymuje całkę na $D_{i}$ używając tylko jednej wartości. Wtedy ${\rm diam}(D_{i})=N^{{-1/d}}$ oraz

$e(f,\overline{M}C_{{d,N}})\,\le\, L\cdot N^{{-(1/2+1/d)}}.$

Ostatecznie, otrzymany w ten sposób wariant losowania warstwowego jest zbieżny z wykładnikiem większym niż $1/2$ . Oczywiście, może to mieć praktyczne znaczenie jedynie dla małych wymiarów $d$ , bowiem dla dużych $d$ wykładnik $1/2+1/d$ jest właściwie równy $1/2$ .

14.3.2. Funkcje kontrolne

Podobny efekt zwiększenia szybkości zbieżności można uzyskać stosując klasyczną Monte Carlo bezpośrednio do funkcji $f-g$ , gdzie $g$ jest pewną specjalnie dobraną funkcją, zwaną funkcją kontrolną.

Rzeczywiście, przedstawmy $f$ w postaci $f=g+(f-g)$ . Wtedy

$I_{d}(f)=I_{d}(g)+I_{d}(f-g),$

co sugeruje zastosowanie następującej metody:

$\widetilde{M}C_{{d,N}}(f)\,:=\, I_{d}(g)+MC_{{d,N}}(f-g).$

Ponieważ

	$\displaystyle I_{d}(f)-\widetilde{M}C_{{d,N}}(f)$	$\displaystyle=$	$\displaystyle(I_{d}(g)+I_{d}(f-g))\,-\,(I_{d}(g)+MC_{{d,N}}(f-g))$
		$\displaystyle=$	$\displaystyle I_{d}(f-g)-MC_{{d,N}}(f-g),$

z twierdzenia 14.1 natychmiast wynika, że błąd

$e(f;\widetilde{M}C_{{d,N}})\,=\, e(f-g;MC_{{d,N}})\,=\,\frac{\sigma(f-g)}{\sqrt{N}}.$

Pozostaje kwestia doboru funkcji $g$ tak, aby istotnie zmniejszyć wariancję $\sigma(f-g)$ . Weźmy najpierw funkcję schodkową postaci

$g(\vec{x})\,=\,\sum _{{i=1}}^{K}c_{i}{\bf 1}_{{D_{i}}}(\vec{x}),$

gdzie

${\bf 1}_{{D_{i}}}(\vec{x})\,=\,\left\{\;\begin{array}[]{cc}1&\quad\vec{x}\in D_{i},\\ 0&\quad\vec{x}\notin D_{i},\end{array}\right.,\qquad 1\le i\le K,$

jest funkcją charakterystyczną zbioru $D_{i}$ , a $c_{i}=f(\vec{x}_{i})$ dla dowolnych $\vec{x}_{i}\in D_{i}$ . Oczywiście, najlepiej byłoby wziąć $\vec{x}_{i}$ tak, aby $c_{i}=I_{d}^{{(i)}}(f)/|D_{i}|$ , ale jest to niemożliwe, bo nie znamy całek $I_{d}^{{(i)}}(f)$ . (Znajomość tych całek nie była konieczna w przypadku losowania warstwowego!) Wtedy dostajemy ten sam efekt jak dla $\overline{M}C_{{d,N}}$ , tzn. dla lipschitzowskiej $f$

$\sigma^{2}(f-g)\,\le\, I_{d}((f-g)^{2})\,\le\, C^{2}\cdot\sum _{{i=1}}^{K}|D_{i}|\,{\rm diam}^{2}(D_{i})$

(14.4)

i dla $g$ odpowiadającej równomiernemu podziałowi na $N$ podkostek otrzymujemy błąd proporcjonalny do $N^{{-(1/2+1/d)}}$ .

W ogólności, funkcję $g$ należy w praktyce wybierać tak, aby dobrze aproksymowała funkcję $f$ . Oczywiście, wybór ten musi bazować na informacji jaką posiadamy o $f$ .

Na przykład, dla funkcji $r$ -gładkich, tzn. dla $f\in{\mathcal{F}}_{r}(D)$ można w ten sposób dostać jeszcze lepszy wykładnik. Rzeczywiście, niech $g=g_{f}$ będzie kawałkami wielomianem stopnia najwyżej $r-1$ po każdej zmiennej interpolującym $f$ , dla równomiernego podziału kostki jednostkowej na $N^{{1/d}}$ podkostek. Z rozdziału 13 wiemy, że wtedy dla funkcji $f\in{\mathcal{F}}_{r}(D)$ błąd interpolacji jest postaci (zob. twierdzenie 13.2)

$|f(\vec{x})-g_{f}(\vec{x})|\,\le\,\frac{C_{{r,d}}B_{r}(f)}{r!}\cdot N^{{-r/d}}$

W konsekwencji, dla tak skonstruowanej metody dostajemy następujący błąd oczekiwany.

Twierdzenie 14.4

Dla $f\in{\mathcal{F}}_{r}(D)$ mamy

$e(f;\widetilde{M}C_{{d,N}})\,\le\,\frac{C_{{r,d}}B_{r}(f)}{r!}\cdot N^{{-(1/2+r/d)}}.$

Dodajmy, że $\widetilde{M}C_{{d,N}}$ jest w istocie metodą mieszaną, gdyż używa wartości funkcji $f$ obliczanych w $N$ punktach wybranych deterministycznie oraz w $N$ punktach wybranych losowo.

Zbieżności $N^{{-(1/2+r/d)}}$ nie da się już poprawić w klasie ${\mathcal{F}}_{r}(D)$ . Dokładniej, można pokazać następujące twierdzenie, które jest odpowiednikiem twierdzenia 13.4 dla algorytmów niedeterministycznych.

Twierdzenie 14.5

Istnieje $c>0$ o następującej własności: dla dowolnej (deterministycznej lub niedeterministycznej) aproksymacji całki wykorzystującej $N$ wartości funkcji istnieje $f\in{\mathcal{F}}_{r}([0,1]^{d})$ dla której $B_{r}(f)=1$ , a błąd oczekiwany aproksymacji całki wynosi co najmniej $c\, N^{{-(1/2+r/d)}}$ .

14.4. Generowanie liczb (pseudo-)losowych

Dotychczas milcząco przyjmowaliśmy, że umiemy generować ciągi

$X_{1},X_{2},X_{3},\ldots,X_{n},\ldots$

niezależnych liczb losowych zgodnie z danym rozkładem prawdopodobieństwa. Nie jest to jednak zadanie trywialne. W praktyce obliczeniowej liczby losowe uzyskujemy przez zastosowanie specjalnych programów. Ponieważ komputer jest urządzeniem deterministycznym, tak uzyskane ciągi nie są idealnie losowe już choćby dlatego, że są okresowe. Fakt ten wpływa na pogorszenie jakości wyniku i w szczególności powoduje, że ich użycie pozwala uzyskać jedynie kilka liczb znaczących, przy czym im większy wymiar $d$ zadania tym gorsza graniczna dokładność. Z tych względów mówimy raczej o generatorach liczb pseudo-losowych.

Generowanie liczb pseudolosowych jest bardzo obszernym tematem, my tylko zwrócimy uwagę na podstawowe metody.

14.4.1. Liniowy generator kongruencyjny

Liniowe generatory kongruencyjne służą generowaniu ciągów losowych $U_{i}$ o rozkładzie jednostajnym na odcinku $[0,1]$ i zdefiniowane są w następujący prosty sposób. Startujemy z $U_{0}=x_{0}$ i kolejno obliczamy dla $i=1,2,3,\ldots$

$\left\{\quad\begin{array}[]{rcl}x_{i}&:=&(a\, x_{{i-1}}+c)\,{\rm mod}\, m,\\ U_{i}&:=&x_{i}/m,\end{array}\right..$

Jakość takiego generatora zależy istotnie o wyboru liczb całkowitych $a$ , $b$ i $m$ . W szczgólności pożądane jest, aby generator miał maksymalny okres $m$ . Jeśli $c\ne 0$ to warunkami dostatecznymi na to są:

(a) $c$ i $m$ są względnie pierwsze,
(b) jeśli $p$ dzieli $m$ to $p$ dzieli $a-1$ ,
(c) jeśli $4$ dzieli $m$ to $4$ dzieli $a-1$ .

Dostęp do dobrego generatora liczb losowych o rozkładnie jednostajnym $U\sim{\rm Unif}([0,1])$ jest ważny również z tego względu, że jest on zwykle podstawą dla konstrukcji generatorów ciągów o bardziej skomplikowanych rozkładach prawdopodobieństwa.

14.4.2. Odwracanie dystrybuanty i ,,akceptuj albo odrzuć”

Jeśli znana jest dystrybuanta żądanego rozkładu, czyli funkcja

$F(x)\,:=\,{\rm Prob}(X\le x),$

oraz łatwo obliczyć jej odwrotność zdefiniowaną jako

$F^{{-1}}(u)\,:=\,\inf\{ x\,:\; F(x)\ge u\},$

to potrzebne ciągi losowe mogą być wygenerowane według wzoru

$X\,:=\, F^{{-1}}(U),\qquad U\sim{\rm Unif}([0,1]).$

Rzeczywiście, mamy

${\rm Prob}(X\le x)\,=\,{\rm Prob}(F^{{-1}}(U)\le x)\,=\,{\rm Prob}(U\le F(x))\,=\, F(x).$

Na przykład, jeśli $F(x)=1-e^{{-x^{2}/2}}$ to można zastosować wzór $X=\sqrt{-2\ln(U)}$ .

Niestety, dla wielu rozkładów dystrybuanta nie może być dokładnie obliczona. Wtedy jakość metody zależy od jakości zastosowanej numerycznej aproksymacji funkcji $F^{{-1}}(x)$ .

Inna uniwersalna metoda generowania liczb losowych o dowolnym rozkładzie na $\mathbb{R}$ , zwana akceptuj albo odrzuć, polega na wykorzystaniu istniejącego ,,dobrego” generatora liczb innego rozkładu na $\mathbb{R}$ . Dokładniej, załóżmy, że dysponujemy generatorem liczb losowych $Y$ zgodnie z rozkładem o gęstości $g$ , a interesują nas liczby $X$ pochodzące z rozkładu o gęstości $f$ . Załóżmy ponadto, że

$f(x)\,\le\, c\cdot g(x)$

dla pewnej stałej $c>0$ . Wtedy możemy użyć następującego algorytmu:

$\{$ .	repeat
.	.	generuj $Y$ zgodnie z $g$ ;
.	.	generuj $U\sim{\rm Unif}([0,1])$
.	until $U\,\le\, f(Y)/(cg(Y))$ ;
.	return $X:=Y$
$\}$

Aby pokazać poprawność takiego generatora zauważmy, że

	$\displaystyle{\rm Prob}(X\in A)$	$\displaystyle=$	$\displaystyle{\rm Prob}(\, Y\in A\,:\; U\le f(Y)/(cg(Y))\;)$
		$\displaystyle=$	$\displaystyle\frac{{\rm Prob}(\, Y\in A,\, U\le f(Y)/(cg(Y))\;)}{P(\, U\le f(Y)/(cg(Y))\;)}.$

Ponieważ dla $a\in[0,1]$ prawdopodobieństwo, że $U\le a$ wynosi $a$ to

${\rm Prob}(U\le f(Y)/(cg(Y))\,)\,=\,\int _{{\mathbb{R}}}\frac{f(x)}{cg(x)}g(x)\, dx\,=\,\frac{1}{c}$

i w konsekwencji dostajemy

	$\displaystyle{\rm Prob}(X\in A)$	$\displaystyle=$	$\displaystyle c\cdot{\rm Prob}(Y\in A,\, U\le f(Y)/(cg(Y))\,)$
		$\displaystyle=$	$\displaystyle c\cdot\int _{A}\frac{f(x)}{cg(x)}g(x)\, dx\,=\,\int _{A}f(x)\, dx.$

14.4.3. Metoda Box-Muller dla rozkładu gaussowskiego

Normalny rozkład gaussowski na $\mathbb{R}$ o funkcji gęstości

$g(t)\,=\,\frac{1}{\sqrt{2\pi}}\, e^{{-t^{2}/2}}$

jest najczęściej stosowanym rozkładem niejednostajnym. Dla rozkładu gaussowskiego bardzo efektywne okazują się algorytmy bazujące na odwracaniu dystrybuanty. Używają one dość skomplikowanych aproksymacji funkcji $F^{{-1}}$ .

Prostszą i najbardziej popularną jest metoda Box-Muller. Generuje ona od razu dwie niezależne liczby $Z_{1}$ i $Z_{2}$ (albo, równoważnie, punkt $(Z_{1},Z_{2})\in\mathbb{R}^{2}$ zgodnie z rozkładem normalnym w $\mathbb{R}^{2}$ ), na podstawie dwóch liczb losowych o rozkładzie jednostajnym.

Przedstawmy $(x,y)\in\mathbb{R}^{2}$ we współrzędnych biegunowych,

$\left\{\begin{array}[]{rcl}x&=&r\cdot\cos\varphi,\\ y&=&r\cdot\sin\varphi,\end{array}\right.$

gdzie $r\in[0,\infty)$ i $\varphi\in[0,2\pi)$ . Metoda polega na wygenerowaniu zmiennych $\varphi$ i $r$ , a następnie zastosowaniu powyższego wzoru. Generowanie $\varphi$ jest proste, bo ma rozkład jednostajny. Policzmy dystrybuantę rokładu zmiennej $r$ . Mamy

	$\displaystyle{\rm Prob}(r\le R)$	$\displaystyle=$	$\displaystyle\frac{1}{2\pi}\int _{{x^{2}+y^{2}\le R^{2}}}e^{{-(x^{2}+y^{2})/2}}\, d(x,y)\,=\,\frac{1}{2\pi}\int _{0}^{{2\pi}}\int _{0}^{R}re^{{-r^{2}/2}}drd\varphi$
		$\displaystyle=$	$\displaystyle\int _{0}^{R}re^{{-r^{2}/2}}dr\,=\,-e^{{-r^{2}/2}}\Big\|_{0}^{R}\,=\, 1-e^{{-R^{2}/2}}.$

Stąd, stosując metodę odwracania dystrybuanty mamy $R=\sqrt{-2\ln U}$ , $U\sim{\rm Unif}([0,1])$ . Rachunki te prowadzą do następującego generatora.

$\{$ .	generuj $U_{1},U_{2}\sim{\rm Unif}([0,1])$ ;
.	$R\,:=\,-2\ln U_{1}$ ; $V\,:=\, 2\pi U_{2}$ ;
.	$Z_{1}\,:=\,\sqrt{R}\,\cos(V)$ ; $Z_{2}\,:=\,\sqrt{R}\,\sin(V)$ ;
.	return $Z_{1},Z_{2}$
$\}$