Zagadnienia

9. Algorytmy Monte Carlo II. Redukcja wariancji

W tym rozdziale omówię niektóre metody redukcji wariancji dla klasycznych algorytmów Monte Carlo, w których losujemy próbki niezależnie, z jednakowego rozkładu. Wiemy, że dla takich algorytmów wariancja estymatora zachowuje się jak const/n. Jedyne, co możemy zrobić – to konstruować takie algorytmy, dla których stała ,,const” jest możliwie mała. Najważniejszą z tych metod faktycznie już poznaliśmy: jest to losowanie istotne, omówione w Rozdziale 8. Odpowiedni wybór ,,rozkładu instrumentalnego” może zmniejszyć wariancję…setki tysięcy razy! Istnieje jeszcze kilka innych, bardzo skutecznych technik. Do podstawowych należą: losowanie warstwowe, metoda zmiennych kontrolnych, metoda zmiennych antytetycznych. Możliwe są niezliczone modyfikacje i kombinacje tych metod. Materiał zawarty w tym rozdziale jest w dużym stopniu zaczerpnięty z monografii Ripleya [18].

9.1. Losowanie warstwowe

Tak jak w poprzednim rozdziale, zedanie polega na obliczeniu wielkości

θ=EπfX=Xfxπxdx, (9.1)

gdzie rozkład prawdopodobieństwa i jego gęstość dla uproszczenia oznaczamy tym samym symbolem π. Losowanie warstwowe polega na tym, że rozbijamy przestrzeń X na sumę k rozłącznych podzbiorów (warstw),

X=h=1kAh,AhAg=hg, (9.2)

i losujemy k niezależnych próbek, po jednej z każdej warstwy. Niech πh będzie gęstostością rozkładu warunkowego zmiennej X przy XAh, czyli

πhx=πxphIxAh,gdzieph=πAh=Ahπxdx. (9.3)

Widać natychmiast, że Bπh(x)dx=π(XB|XAh). Rozbijamy teraz całkę, którą chcemy obliczyć:

θ=hphfxπhxdx. (9.4)

Możemy użyć następującego estymatora warstwowego:

θnstra=hphnhi=1nhfXhi, (9.5)

gdzie

Xh1,,Xhnhi.i.dπh, (9.6)

jest próbką rozmiaru nh wylosowaną z h-tej warstwy (h=1,,k). Jest to estymator nieobciążony,

Varθnstra=hph2nhσh2, (9.7)

gdzie σh2=Varπ(f(X)|xAh).

Porównajmy estymator warstwowy (9.5) ze ,,zwykłym” estymatorem

θnCMC=1ni=1nfXi. (9.8)

opartym na jednej próbce

X1,,Xni.i.dπ. (9.9)

Oczywiście

VarθnCMC=1nσ2, (9.10)

gdzie σ2=VarπfX. Żeby porównanie było ,,uczciwe” rozważmy próbkę liczności n=hnh. Jeśli decydujemy się na sumaryczną liczność próbki n, to możemy w różny sposób ,,rozdzielić” to n pomiędzy warstwami. To się nazywa alokacja próbki.

Stwierdzenie 9.1 (Alokacja proporcjonalna)

Jeżeli nh=phn dla n=1,,k, to

Varθnstra=1nhphσh21nσ2=VarθnCMC. (9.11)
Dowód

Wyrażenie na wariancję wynika znatychmiast z podstawienia nh=phn w ogólnym wzorze (9.7). Nierówność wynika z następującej tożsamości:

σ2=hphσh2+hphθh-θ2, (9.12)

gdzie θh=f(x)πh(x)dx=Eπ(f(X)|XAh). Jest to klasyczny wzór na dekompozycję wariancji na ,,wariancję wewnątrz warstw” (pierwszy składnik w (9.12)) i ,,wariancję pomiędzy warstwami” (drugi składnik w (9.12)). Zdefiniujemy zmienną losową H jako ,,numer warstwy do której wpada Xπ”, czyli {H=h}={XAh}. Możemy teraz wzór (9.12) przepisać w dobrze znanej postaci

Varf(X)=EVar(f(X)|H)+VarE(f(X)|H). (9.13)

Wzór (9.12) podpowiada, jak dzielić przestrzeń na warstwy. Największy zysk w porównaniu z losowaniem nie-warstwowym jest wtedy, gdy ,,wariancja międzywarstwowa” jest dużo większa od ,,wewnątrzwarstwowej” Warstwy należy więc wybierać tak, żeby funkcja πcfx była możliwie bliska stałej na każdym zbiorze Ah i różniła się jak najbardziej pomiędzy poszczególnymi zbiorami.

Stwierdzenie 9.1 pokazuje, że zawsze zyskujemy na losowaniu warstwowym, jeśli zastosujemy najprostszą, proporcjonalną alokację. Okazuje się, że nie jest to alokacja najlepsza. Jerzy Spława-Neyman odkrył prostą regułę wyznaczania alokacji optymalnej. Wychodzimy od wzoru (9.7) i staramy się zoptymalizować prawą stronę przy warunku n=hnh. Poszukujemy zatem rozwiązania zadania

hph2nhσh2=min!hnh-n=0 (9.14)

(względem zmiennych nh). Zastosujmy metodę mnożników Lagrange'a. Szukamy minimum

L=hph2nhσh2+λhnh-n=min! (9.15)

Obliczamy pochodną i przyrównujemy do zera:

nhL=-ph2nh2σh2+λ=0. (9.16)

Stąd natychmiast otrzymujemy rozwiązanie: nhσhph.

Stwierdzenie 9.2 (Alokacja optymalna, J. Neyman)

Estymator warstwowy ma najmniejszą wariancję jeśli alokacja n losowanych punktów jest dana wzorem

nh=σhphgσgpg,h=1,,k. (9.17)

Zignorowaliśmy tutaj niewygodne wymaganie, że liczności próbek nh muszą być całkowite. Gdyby to wziąć pod uwagę, rozwiązanie stałoby się skomplikowane, a zysk praktyczny z tego byłby znikomy. W praktyce rozwiązanie neymanowskie zaokrągla się do liczb całkowitych i koniec. Ważniejszy jest inny problem. Żeby wyznaczyć alokację optymalną, trzeba znać nie tylko prawdopodobieństwa ph ale i wariancje warstwowe σh2. W praktyce często opłaca się wylosować wstępne próbki, na podstawie których estymuje się wariancje σh2. Dopiero potem alokuje się dużą, roboczą próbkę rozmiaru n, która jest wykorzystana do obliczania docelowej całki.

9.2. Zmienne kontrolne

Idea zmiennych kontrolnych polega na rozbiciu docelowej całki (wartości oczekiwanej) na dwa składniki, z których jeden umiemy obliczyć analitycznie. Metodę Monte Carlo stosujemy do drugiego składnika. Przedstawmy wielkość obliczaną w postaci

θ=EπfX=Xfxπxdx=Xfx-kxπxdx+Xkxπxdx. (9.18)

Dążymy do tego, żeby całka funkcji k była obliczona analitycznie (lub numerycznie) a różnica f-k była możliwie bliska stałej, bo wtedy wariancja metody Monte Carlo jest mała. Funkcję k lub zmienną losową kX nazywamy zmienną kontrolną.

Dla uproszczenia połóżmy Y=fX, gdzie Xπ. Przypuśćmy, że zmiennej kontrolnej będziemy szukać pośród kombinacji liniowych funkcji k1,,kd o znanych całkach. Niech Zj=kjX i Z=Z1,,Zd. Przy tym stale pamiętajmy, że Xπ i będziemy pomijać indeks π przy wartościach oczekiwanych i wariancjach. Zatem

kx=j=1dβjkjx. (9.19)

Innymi słowy poszukujemy wektora współczynników β=β1,,βd który minimalizuje wariancję VarY-βZ. Wykorzystamy następujący standardowy wynik z teorii regresji liniowej.

Stwierdzenie 9.3

Niech Y i Z będą zmiennymi losowymi o skończonych drugich momentach, wymiaru odpowiednio 1 i d. Zakładamy dodatkowo odwracalność macierzy wariancji-kowariancji VARZ. Kowariancję pomiędzy Y i Z traktujemy jako wektor wierszowy i oznaczamy przez COVY,Z Spośród zmiennych losowych postaci Y-βZ, najmniejszą wariancję otrzymujemy dla

β*=COVY,ZVARZ-1. (9.20)

Ta najmniejsza wartość wariancji jest równa

VarY-β*Z=VarY-COVY,ZVARZ-1COVZ,Y. (9.21)

Przepiszmy tem wynik w bardziej sugestywnej formie. Niech VarY=σ2. Można pokazać, że β* maksymalizuje korelację pomiędzy Y i βZ. Napiszmy

ϱY.Z=maxβcorrY,βZ=corrY,β*Z=COVY,ZVARZ-1COVZ,YVarY. (9.22)

Niech teraz θncontr będzie estymatorem zmiennych kontrolnych. To znaczy, że losujemy próbkę X1,,Xnπ,

θncontr=1nYi-β*Zi+β*μ, (9.23)

gdzie Zi=Zi1,,Zin=f1Xi,,fdXi, zaś μj=EfjX są obliczone analitycznie lub numerycznie μ=μ1,,μd. Wariancja estymatora jest wyrażona wzorem

Varθncontr=1nσ21-ϱY.Z2, (9.24)

co należy porównać z wariancją σ2/n ,,zwykłego” estymatora.

Zróbmy podobną uwagę, jak w poprzednim podrozdziale. Optymalny wybór współczynników regresji wymaga znajomości wariancji i kowariancji, których obliczenie może być …(i jest zazwyczaj!) trudniejsze niż wyjściowe zadanie obliczenia wartości oczekiwanej. Niemniej, można najpierw wylosować wstępną próbkę, wyestymować potrzebne wariancje i kowariancje (nawet niezbyt dokładnie) po to, żeby dla dużej, roboczej próbki skonstruować dobre zmienne kontrolne.

Wspomnijmy na koniec, że dobieranie zmiennej kontrolnej metodą regresji liniowej nie jest jedynym sposobem. W konkretnych przykładach można spotkać najróżniejsze, bardzo pomysłowe konstrukcje, realizujące podstawową ideę zmiennych kontrolnych.

9.3. Zmienne antytetyczne

Przypuśćmy, że estymujemy wielkość θ=EπfX. Jeśli mamy dwie zmienne losowe X i X o jednakowym rozkładzie π ale nie zakładamy ich niezależności, to

VarfX+fX2=12VarX1+corrfX,fX=12σ21+ϱ.

Jeśli ϱ=corrfX,fX<0, to wariancja w powyższym wzorze jest mniejsza niż σ2/2, czyli mniejsza niż w przypadku niezależności X i X. To sugeruje, żeby zamiast losować niezależnie n zmiennych X1,,Xn, wygenerować pary zmiennych ujemnie skorelowanych. Załóżmy, że n=2k i mamy k niezależnych par X1,X1,,Xk,Xk, a więc łącznie n zmiennych. Porównajmy wariancję ,,zwykłego” estymatora θnCMC i estymatora θnant wykorzystującego ujemne skorelowanie par:

θnCMC=1ni=1nfXi,VarθnCMC=σ2nθnant=1ni=1n/2fXi=fXi,Varθnant=σ2n1+ϱ.

Wariancja estymatora θnant jest tym mniejsza, im ϱ bliższe -1 (im bardziej ujemnie skorelowane są pary). Standardowym sposobem generowania ujemnie skorelowanych par jest odwracanie dystrybuanty z użyciem ,,odwróconych loczb losowych”.

Stwierdzenie 9.4

Jeśli h:]0,1[R jest funkcją monotoniczną różną od stałej, 01hu2du< i UU0,1, to

CovhU,h1-U<0. (9.25)
Dowód

Bez straty ogólności załóżmy, że h jest niemalejąca. Niech μ=EhU=01hudu i t=supu:h1-u>μ. Łatwo zauważyć, że 0<t<1. Zauważmy, że

CovhU,h1-U=EhUh1-U-μ=01huh1-u-μdu=0thuh1-u-μdu+t1huh1-u-μdu<0thth1-u-μdu+t1hth1-u-μdu=ht01h1-u-μdu=0, (9.26)

ponieważ dla 0<u<t mamy h1-u-μ>0 i dla t<u<1 mamy h1-u-μ0.

Przypomnijmy, że dla dowolnej dystrybuanty G określamy uogólnioną funkcję odwrotną G- następującym wzorem:

G-u=infx:Gxu. (9.27)

Natychmiast wnioskujemy ze Stwierdzenia 9.4, że CovG-U,G-1-U<0. W ten sposób możemy produkować ujemnie skorelowane pary zmiennych o zadanej dystrybuancie. Okazuje się, że są to najbardziej ujemnie skorelowane pary. Jeśli XG i XG, to

CovX,XCovG-U,G-1-U. (9.28)

Powyższy fakt ma oczywiste znaczenie z punktu widzenia metod Monte Carlo. Udowodnimy ogólniejsze twierdzenie.

Twierdzenie 9.1

Jeżeli XF, YG oraz EX2<, EY2< to

CovF-U,G-1-UCovX,YCovF-U,G-U. (9.29)

Twierdzenie 9.2 wynika z trzech poniższych faktów. Każdy z nich jest sam w sobie interesujący. Zaczniemy od sławnego wyniku Frecheta.

Twierdzenie 9.2 (Ograniczenia Frecheta)

Jeżeli PXx=Fx, PYx=Gy oraz PXx,Yy=Hx,y oznaczają, odpowiednio, dystrybuanty brzegowe oraz łączną dystrybuantę dwóch zmiennych losowych to

max0,Fx+Gy-1Hx,yminFx,Gy. (9.30)

Istnieją rozkłady łączne o brzegowych F i G, dla których jest osiągane ograniczenie dolne i ograniczenie górne.

Dowód

Ograniczenie górne wynika z oczywistych nierówności

P(Xx,Yy)P(Xx)=F(x),P(Xx,Yy)P(Yy)=G(y). (9.31)

Ograniczenie dolne jest równie proste:

P(Xx,Yy)=P(Xx)-P(Xx,Y>y)P(Xx)-P(Y>y)=F(x)-[1-G(y)]. (9.32)

Pozostała do pokazania osiągalność. Następujący lemat jest jednym z piękniejszych przykładów ,,symulacyjnego punktu widzenia” w rachunku prawdopodobieństwa.

Lemat 9.1

Jeśli UU0,1 to

P(F-(U)x,G-(U)y)=min(F(x),G(y));P(F-(U)x,G-(1-U)y)=max(0,F(x)+G(y)-1). (9.33)
Dowód

Pierwsza równość jest oczywista:

P(F-(U)x,G-(U)y)=P(UF(x),UG(y))=min(F(x),G(y)). (9.34)

Druga równość też jest oczywista:

P(F-(U)x,G-(1-U)y)=P(UF(x),1-UG(y))=P(1-G(y)UF(x))=max(0,F(x)-[1-G(y)]). (9.35)

Głębokie twierdzenie Frecheta składa się więc z kilku dość oczywistych spostrzeżeń. Zeby udowodnić Twierdzenie 9.1 potrzeba jeszcze jednego ciekawego lematu.

Lemat 9.2

Niech Fx, Gy i Hx,y oznaczają, odpowiednio, dystrybuanty brzegowe oraz łączną dystrybuantę zmiennych losowych X i Y. Jeśli EX2<, EY2< to

CovX,Y=Hx,y-FxFydxdy. (9.36)
Dowód

Niech X1,Y1 i X2,Y2 będą niezależnymi parami o jednakowym rozkładzie takim jak para X,Y. Wtedy

2Cov(X,Y)=E(X1-X2)(Y1-Y2)=E[I(X1x)-I(X2x)][I(Y1y)-I(Y2y)]dxdy=E[I(X1x)-I(X2x)][I(Y1y)-I(Y2y)]dxdy=[P(X1x,Y1y)+P(X2x,Y2y)-P(X1x,Y2y)-P(X2x,Y1y)]dxdy=[2P(Xx,Yy)-2P(Xx)P(Yy)]dxdy. (9.37)

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.