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

Id⁢f⁢:=⁢∫0,1df⁢x→⁢⁢d⁢x→⁢=⁢∫01∫01⋯⁢∫01︸d⁢f⁢x1,x2,…,xd⁢⁢d⁢x1⁢d⁢x2⁢⋯⁢d⁢xd.

Zakładamy przy tym, że f:0,1d→R jest funkcją, której kwadrat jest całkowalny,

∫0,1df⁢x→2⁢⁢d⁢x→⁢<⁢∞.

Konsekwencją zastosowania losowości jest to, że przy różnych realizacjach metody otrzymujemy różne wyniki, w zależności od wyboru punktów t→j. Wynik M⁢Cd,N⁢f jest więc zmienną losową, której wartość oczekiwana wynosi

	E⁢M⁢Cd,N⁢f	=	∫0,1d⋅NM⁢Cd,N⁢f;t→1,…,t→N⁢⁢d⁢t→1⁢⋯⁢d⁢t→N
		=	1N⁢∑j=1N∫0,1df⁢t→⁢⁢d⁢t→⁢=⁢Id⁢f.

Ponieważ różnica Id⁢f-M⁢Cd,N⁢f jest też zmienną losową, za błąd metody Monte Carlo dla danej funkcji f przyjmiemy odchylenie standardowe,

e⁢f;M⁢Cd,N⁢:=⁢E⁢Id⁢f-M⁢Cd,N⁢f2.

Zanim przystąpimy do dowodu zauważmy, że σ⁢f jest dobrze zdefiniowaną wielkością, bowiem nierówność

Id⁢f≤Id⁢f2

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

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 c1<c2 mamy

limN→∞Prob(c1⁢σ⁢fN≤Id(f)-MCd,N(f)≤c2⁢σ⁢fN)=12⁢π∫c1c2e-t2/2dt,

gdzie Prob oznacza prawdopodobieństwo względem rozkładu jednostajnego na 0,1d⁢N.

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

Id,ω(f):=∫Rdf(x→)ω(x→)dx→,gdzie∫Rdω(x→)dx→=:W<∞,

możemy bowiem zastosować wzór

M⁢Cd,Nω⁢f⁢:=⁢WN⋅∑j=1Nf⁢t→j,

przy czym t→1,…,t→N są tym razem punktami wybranymi losowo i niezależnie od siebie, zgodnie z rozkładem na Rd o gęstości ω/W.

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

Podzielmy obszar całkowania D na K rozłącznych podzbiorów Di tak, że

D=⋃i=1KDi

i zastosujmy Monte Carlo do całkowania po każdym Di, tzn. całkę ∫Df⁢x→⁢⁢d⁢x→ przybliżymy wielkością

M¯⁢Cd,N⁢f⁢:=⁢∑i=1KM⁢Cd,Nii⁢f,

gdzie M⁢Cd,Nii jest metodą Monte Carlo zastosowaną do całki

Idi⁢f⁢:=⁢∫Dif⁢x→⁢⁢d⁢x→,⁢1≤i≤K,

oraz N=∑i=1KNi.

Oznaczmy przez Di objętość d-wymiarową podzbioru Di. Ponieważ zmienne losowe Idi⁢f-M⁢Cd,Nii⁢f są parami niezależne dla 1≤i≤K, na podstawie twierdzenia 14.2 mamy

	E⁢Id⁢f-M¯⁢Cd,N⁢f2	=	E⁢∑i=1KIdi⁢f-M⁢Cd,Ni⁢f2
		=	∑i=1KE⁢Idi⁢f-M⁢Cd,Nii⁢f2
		=	∑i=1K1Ni⁢Di⁢⁢Idi⁢f2-Idi⁢f2.

Przyjmijmy teraz, że

Ni⁢=⁢Di⋅N,⁢1≤i≤K,

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

e⁢f,M¯⁢Cd,N⁢=⁢1N⋅Id⁢f2-∑i=1K1Di⁢Idi⁢f2.

(14.2)

Błąd tak zdefiniowanej metody M¯⁢Cd,N nie jest większy od błędu klasycznej metody M⁢Cd,N z Twierdzenia 14.1.

Widzimy, że stosując losowanie warstwowe z ustalonym podziałem na K podzbiory Di 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 M¯⁢Cd,N do postaci

e⁢f,M¯⁢Cd,N⁢=⁢1N⁢∑i=1KIdi⁢f-ci2.

(14.3)

gdzie

ci:=Idi⁢fDi,⁢1≤i≤K.

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

f⁢x→-f⁢y→⁢≤⁢L⋅x→-y→∞,⁢x→,y→∈D.

Wtedy istnieją x→i∈D¯i takie, że f⁢x→i=ci, a stąd i z lipschitzowskości f mamy, że dla dowolnego x→∈Di

f⁢x→-ci⁢≤⁢L⋅x→-x→i∞⁢≤⁢L⋅diam∞⁢Di,

gdzie

diam∞⁢Di⁢:=⁢sup⁡x→-y→∞⁢:⁢x→,y→∈Di

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

e⁢f,M¯⁢Cd,N⁢≤⁢LN⁢∑i=1KDi⁢⁢diam2⁢Di.

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

e⁢f,M¯⁢Cd,N⁢≤⁢L⋅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.

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

Id⁢f=Id⁢g+Id⁢f-g,

co sugeruje zastosowanie następującej metody:

M~⁢Cd,N⁢f⁢:=⁢Id⁢g+M⁢Cd,N⁢f-g.

Ponieważ

	Id⁢f-M~⁢Cd,N⁢f	=	Id⁢g+Id⁢f-g⁢-⁢Id⁢g+M⁢Cd,N⁢f-g
		=	Id⁢f-g-M⁢Cd,N⁢f-g,

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

e⁢f;M~⁢Cd,N⁢=⁢e⁢f-g;M⁢Cd,N⁢=⁢σ⁢f-gN.

Pozostaje kwestia doboru funkcji g tak, aby istotnie zmniejszyć wariancję σ⁢f-g. Weźmy najpierw funkcję schodkową postaci

g⁢x→⁢=⁢∑i=1Kci⁢1Di⁢x→,

gdzie

1Di(x→)={1⁢x→∈Di,0⁢x→∉Di,,1≤i≤K,

jest funkcją charakterystyczną zbioru Di, a ci=f⁢x→i dla dowolnych x→i∈Di. Oczywiście, najlepiej byłoby wziąć x→i tak, aby ci=Idi⁢f/Di, ale jest to niemożliwe, bo nie znamy całek Idi⁢f. (Znajomość tych całek nie była konieczna w przypadku losowania warstwowego!) Wtedy dostajemy ten sam efekt jak dla M¯⁢Cd,N, tzn. dla lipschitzowskiej f

σ2⁢f-g⁢≤⁢Id⁢f-g2⁢≤⁢C2⋅∑i=1KDi⁢⁢diam2⁢Di

(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∈Fr⁢D można w ten sposób dostać jeszcze lepszy wykładnik. Rzeczywiście, niech g=gf 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 N1/d podkostek. Z rozdziału 13 wiemy, że wtedy dla funkcji f∈Fr⁢D błąd interpolacji jest postaci (zob. twierdzenie 13.2)

f⁢x→-gf⁢x→⁢≤⁢Cr,d⁢Br⁢fr!⋅N-r/d

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

Dodajmy, że M~⁢Cd,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 Fr⁢D. Dokładniej, można pokazać następujące twierdzenie, które jest odpowiednikiem twierdzenia 13.4 dla algorytmów niedeterministycznych.

Liniowe generatory kongruencyjne służą generowaniu ciągów losowych Ui o rozkładzie jednostajnym na odcinku 0,1 i zdefiniowane są w następujący prosty sposób. Startujemy z U0=x0 i kolejno obliczamy dla i=1,2,3,…

{xi:=a⁢⁢xi-1+c⁢⁢mod⁢⁢m,Ui:=xi/m,.

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≠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∼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.

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

F⁢x⁢:=⁢ProbX≤x,

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

F-1⁢u⁢:=⁢inf⁡x⁢:⁢F⁢x≥u,

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

X⁢:=⁢F-1⁢U,⁢U∼Unif⁢0,1.

Rzeczywiście, mamy

Prob(X≤x)=Prob(F-1(U)≤x)=Prob(U≤F(x))=F(x).

Na przykład, jeśli F⁢x=1-e-x2/2 to można zastosować wzór X=-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 R, zwana akceptuj albo odrzuć, polega na wykorzystaniu istniejącego ,,dobrego” generatora liczb innego rozkładu na 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⁢≤⁢c⋅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∼Unif⁢0,1
.	until U⁢≤⁢f⁢Y/c⁢g⁢Y;
.	return X:=Y
}

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

	ProbX∈A	=	Prob(Y∈A:U≤f(Y)/(cg(Y)))
		=	Prob⁢Y∈A,⁢U≤f⁢Y/c⁢g⁢Y⁢P⁢U≤f⁢Y/c⁢g⁢Y⁢.

Ponieważ dla a∈0,1 prawdopodobieństwo, że U≤a wynosi a to

Prob(U≤f(Y)/(cg(Y)))=∫Rf⁢xc⁢g⁢xg(x)dx=1c

i w konsekwencji dostajemy

	ProbX∈A	=	c⋅ProbY∈A,⁢U≤f⁢Y/c⁢g⁢Y⁢
		=	c⋅∫Af⁢xc⁢g⁢x⁢g⁢x⁢⁢d⁢x⁢=⁢∫Af⁢x⁢⁢d⁢x.

Normalny rozkład gaussowski na R o funkcji gęstości

g⁢t⁢=⁢12⁢π⁢⁢e-t2/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 Z1 i Z2 (albo, równoważnie, punkt Z1,Z2∈R2 zgodnie z rozkładem normalnym w R2), na podstawie dwóch liczb losowych o rozkładzie jednostajnym.

Przedstawmy x,y∈R2 we współrzędnych biegunowych,

⁡x=r⋅cos⁡φ,y=r⋅sin⁡φ,⁢

gdzie r∈0,∞ i φ∈0,2⁢π. Metoda polega na wygenerowaniu zmiennych φ i r, a następnie zastosowaniu powyższego wzoru. Generowanie φ jest proste, bo ma rozkład jednostajny. Policzmy dystrybuantę rokładu zmiennej r. Mamy

	Probr≤R	=	12⁢π⁢∫x2+y2≤R2e-x2+y2/2⁢⁢d⁢x,y⁢=⁢12⁢π⁢∫02⁢π∫0Rr⁢e-r2/2⁢d⁢r⁢d⁢φ
		=	∫0Rr⁢e-r2/2⁢d⁢r⁢=-e-r2/20R⁢= 1-e-R2/2.

Stąd, stosując metodę odwracania dystrybuanty mamy R=-2⁢ln⁡U, U∼Unif⁢0,1. Rachunki te prowadzą do następującego generatora.

{ .	generuj U1,U2∼Unif⁢0,1;
.	R⁢:=-2⁢ln⁡U1;  V⁢:= 2⁢π⁢U2;
.	Z1⁢:=⁢R⁢⁢cos⁡V;  Z2⁢:=⁢R⁢⁢sin⁡V;
.	return Z1,Z2
}