Zagadnienia

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ę

Idf:=0,1dfxdx=010101dfx1,x2,,xddx1dx2dxd.

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

0,1dfx2dx<.
Definicja 14.1

Klasyczna metoda Monte Carlo polega na przybliżeniu Idf średnią arytmetyczną wartości funkcji f w losowo wybranych punktach, tzn.

MCd,Nf=MCd,Nf;t1,t2,,tN:=1Nj=1Nftj,

gdzie t1,t2,,tN są punktami wylosowanymi niezależnie od siebie, każdy zgodnie z rozkładem jednostajnym na 0,1d.

Konsekwencją zastosowania losowości jest to, że przy różnych realizacjach metody otrzymujemy różne wyniki, w zależności od wyboru punktów tj. Wynik MCd,Nf jest więc zmienną losową, której wartość oczekiwana wynosi

EMCd,Nf=0,1dNMCd,Nf;t1,,tNdt1dtN
=1Nj=1N0,1dftdt=Idf.

Ponieważ różnica Idf-MCd,Nf jest też zmienną losową, za błąd metody Monte Carlo dla danej funkcji f przyjmiemy odchylenie standardowe,

ef;MCd,N:=EIdf-MCd,Nf2.
Twierdzenie 14.1

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

ef;MCd,N=σfN,

gdzie

σf:=Idf2-Id2f

jest wariancją funkcji f.

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

IdfIdf2

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

Dowód

Oznaczmy, dla uproszczenia, zmienną losową X=MCd,Nf. Wtedy

EX-EX2=EXX-EX-EXX-EX=EX2-E2X.(14.1)

Ponadto

EX2=E1Nj=1Nftj2=1N2Ej=1Nf2tj+ijftiftj
=1N2NIdf2+N2-NId2f=1NIdf2+1-1NId2f,

gdzie skorzystaliśmy z niezależności zmiennych losowych fti i ftj dla ij. Stąd i z (14.1) dostajemy

e2f;MCd,N=EX-EX2=1NIdf2+1-1NId2f-Id2f=1NIdf2-Id2f,

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

limNProb(c1σfNId(f)-MCd,N(f)c2σfN)=12πc1c2e-t2/2dt,

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

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

Id,ω(f):=Rdf(x)ω(x)dx,gdzieRdω(x)dx=:W<,

możemy bowiem zastosować wzór

MCd,Nωf:=WNj=1Nftj,

przy czym t1,,tN 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.

Twierdzenie 14.2

Niech Id,ωf2=Rdf2xωxdx<. Wtedy

ef;MCd,Nω=σωfN,

gdzie

σωf=WId,ωf2-Id,ω2f.

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 σωfN-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 σω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,1d.

14.3.1. Losowanie warstwowe

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ę Dfxdx przybliżymy wielkością

M¯Cd,Nf:=i=1KMCd,Niif,

gdzie MCd,Nii jest metodą Monte Carlo zastosowaną do całki

Idif:=Difxdx,1iK,

oraz N=i=1KNi.

Oznaczmy przez Di objętość d-wymiarową podzbioru Di. Ponieważ zmienne losowe Idif-MCd,Niif są parami niezależne dla 1iK, na podstawie twierdzenia 14.2 mamy

EIdf-M¯Cd,Nf2=Ei=1KIdif-MCd,Nif2
=i=1KEIdif-MCd,Niif2
=i=1K1NiDiIdif2-Idif2.

Przyjmijmy teraz, że

Ni=DiN,1iK,

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

ef,M¯Cd,N=1NIdf2-i=1K1DiIdif2.(14.2)

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

Twierdzenie 14.3

Dla dowolnej funkcji f takiej, że Idf2< mamy

ef,M¯Cd,Nef,MCd,N,

przy czym równość zachodzi tylko wtedy gdy iloraz Idif/Di jest stały, niezależnie od i.

Dowód

Rzeczywiście, oznaczając

ai=Di,bi=IdifDi,

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

Id2f=i=1Kaibi2i=1Kai2i=1Kbi2=i=1Kbi2=i=1K1DiIdif2,

przy czym równość zachodzi tylko wtedy gdy wektory a1,,aKT i b1,,bKT 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 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

ef,M¯Cd,N=1Ni=1KIdif-ci2.(14.3)

gdzie

ci:=IdifDi,1iK.

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

fx-fyLx-y,x,yD.

Wtedy istnieją xiD¯i takie, że fxi=ci, a stąd i z lipschitzowskości f mamy, że dla dowolnego xDi

fx-ciLx-xiLdiamDi,

gdzie

diamDi:=supx-y:x,yDi

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

ef,M¯Cd,NLNi=1KDidiam2Di.

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 diamDi=N-1/d oraz

ef,M¯Cd,NLN-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

Idf=Idg+Idf-g,

co sugeruje zastosowanie następującej metody:

M~Cd,Nf:=Idg+MCd,Nf-g.

Ponieważ

Idf-M~Cd,Nf=Idg+Idf-g-Idg+MCd,Nf-g
=Idf-g-MCd,Nf-g,

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

ef;M~Cd,N=ef-g;MCd,N=σf-gN.

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

gx=i=1Kci1Dix,

gdzie

1Di(x)={1xDi,0xDi,,1iK,

jest funkcją charakterystyczną zbioru Di, a ci=fxi dla dowolnych xiDi. Oczywiście, najlepiej byłoby wziąć xi tak, aby ci=Idif/Di, ale jest to niemożliwe, bo nie znamy całek Idif. (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

σ2f-gIdf-g2C2i=1KDidiam2Di(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 fFrD 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 fFrD błąd interpolacji jest postaci (zob. twierdzenie 13.2)

fx-gfxCr,dBrfr!N-r/d

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

Twierdzenie 14.4

Dla fFrD mamy

ef;M~Cd,NCr,dBrfr!N-1/2+r/d.

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 FrD. 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 fFr0,1d dla której Brf=1, a błąd oczekiwany aproksymacji całki wynosi co najmniej cN-1/2+r/d.

14.4. Generowanie liczb (pseudo-)losowych

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

X1,X2,X3,,Xn,

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 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:=axi-1+cmodm,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 c0 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 UUnif0,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

Fx:=ProbXx,

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

F-1u:=infx:Fxu,

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

X:=F-1U,UUnif0,1.

Rzeczywiście, mamy

Prob(Xx)=Prob(F-1(U)x)=Prob(UF(x))=F(x).

Na przykład, jeśli Fx=1-e-x2/2 to można zastosować wzór X=-2lnU.

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-1x.

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

fxcgx

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

{ . repeat
. . generuj Y zgodnie z g;
. . generuj UUnif0,1
. until UfY/cgY;
. return X:=Y
}

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

ProbXA=Prob(YA:Uf(Y)/(cg(Y)))
=ProbYA,UfY/cgYPUfY/cgY.

Ponieważ dla a0,1 prawdopodobieństwo, że Ua wynosi a to

Prob(Uf(Y)/(cg(Y)))=Rfxcgxg(x)dx=1c

i w konsekwencji dostajemy

ProbXA=cProbYA,UfY/cgY
=cAfxcgxgxdx=Afxdx.

14.4.3. Metoda Box-Muller dla rozkładu gaussowskiego

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

gt=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,Z2R2 zgodnie z rozkładem normalnym w R2), na podstawie dwóch liczb losowych o rozkładzie jednostajnym.

Przedstawmy x,yR2 we współrzędnych biegunowych,

x=rcosφ,y=rsinφ,

gdzie r0, 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

ProbrR=12πx2+y2R2e-x2+y2/2dx,y=12π02π0Rre-r2/2drdφ
=0Rre-r2/2dr=-e-r2/20R= 1-e-R2/2.

Stąd, stosując metodę odwracania dystrybuanty mamy R=-2lnU, UUnif0,1. Rachunki te prowadzą do następującego generatora.

{ . generuj U1,U2Unif0,1;
. R:=-2lnU1;  V:= 2πU2;
. Z1:=RcosV;  Z2:=RsinV;
. return Z1,Z2
}

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.