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

θ=Eπ⁢f⁢X=∫Xf⁢x⁢π⁢x⁢d⁢x,

(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,⁢Ah∩Ag=∅⁢h≠g,

(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 X∈Ah, czyli

πh⁢x=π⁢xph⁢Ix∈Ah,⁢gdzie⁢⁢ph=π⁢Ah=∫Ahπ⁢x⁢d⁢x.

(9.3)

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

θ=∑hph⁢∫f⁢x⁢πh⁢x⁢d⁢x.

(9.4)

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

θ⌃nstra=∑hphnh⁢∑i=1nhf⁢Xh⁢i,

(9.5)

gdzie

Xh⁢1,…,Xh⁢nh∼i.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)|x∈Ah).

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

θ⌃nCMC=1n⁢∑i=1nf⁢Xi.

(9.8)

opartym na jednej próbce

X1,…,Xn∼i.i.dπ.

(9.9)

Oczywiście

Var⁢θ⌃nCMC=1n⁢σ2,

(9.10)

gdzie σ2=Varπ⁢f⁢X. Ż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.

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 π⁢c⁢f⁢x 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:

∂∂⁡nh⁢L=-ph2nh2⁢σh2+λ=0.

(9.16)

Stąd natychmiast otrzymujemy rozwiązanie: nh∝σh⁢ph.

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.

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π⁢f⁢X=∫Xf⁢x⁢π⁢x⁢d⁢x=∫Xf⁢x-k⁢x⁢π⁢x⁢d⁢x+∫Xk⁢x⁢π⁢x⁢d⁢x.

(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ą k⁢X nazywamy zmienną kontrolną.

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

k⁢x=∑j=1dβj⁢kj⁢x.

(9.19)

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

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

ϱY.Z=maxβ⁡corr⁢Y,β⊤⁢Z=corr⁢Y,β*⊤⁢Z=COV⁢Y,Z⁢VAR⁢Z-1⁢COV⁢Z,YVar⁢Y.

(9.22)

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

θ⌃ncontr=1n⁢∑Yi-β*⊤⁢Zi+β*⊤⁢μ,

(9.23)

gdzie Zi⊤=Zi⁢1,…,Zi⁢n=f1⁢Xi,…,fd⁢Xi, zaś μj=E⁢fj⁢X są obliczone analitycznie lub numerycznie μ=μ1,…,μd⊤. Wariancja estymatora jest wyrażona wzorem

Var⁢θ⌃ncontr=1n⁢σ2⁢1-ϱ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.

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

Var⁢f⁢X+f⁢X′2=12⁢Var⁢X⁢1+corr⁢f⁢X,f⁢X′=12⁢σ2⁢1+ϱ.

Jeśli ϱ=corr⁢f⁢X,f⁢X′<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=2⁢k 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⁢=1n⁢∑i=1nf⁢Xi,⁢Var⁢θ⌃nCMC=σ2n⁢⁢θ⌃nant⁢=1n⁢∑i=1n/2f⁢Xi=f⁢Xi′,⁢Var⁢θ⌃nant=σ2n⁢1+ϱ.⁢

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”.

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

G-⁢u=inf⁡x:G⁢x≥u.

(9.27)

Natychmiast wnioskujemy ze Stwierdzenia 9.4, że Cov⁢G-⁢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 X∼G i X′∼G, to

Cov⁢X,X′≥Cov⁢G-⁢U,G-⁢1-U.

(9.28)

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

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.

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.

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.