Zagadnienia

5.1 Ogólne metody
- 5.1.1 Metoda rozkładów warunkowych
- 5.1.2 Metoda przeksztalceń
5.2 Kilka ważnych przykładów
- 5.2.1 Rozkłady sferyczne i eliptyczne
- 5.2.2 Rozklady Dirichleta

5. Generowanie zmiennych losowych III. Rozkłady wielowymiarowe

5.1. Ogólne metody

Wiele spośród ogólnych metod generowania zmiennych losowych jest całkowicie niezależnych od wymiaru. W szczególności, metody eliminacji, kompozycji i przekształceń są z powodzeniem stosowane do generowania zmiennych losowych wielowymiarowych. Wyjątek stanowi ,,najbardziej ogólna” metoda odwracania dystrybuanty, która nie ma naturalnego odpowiednika dla wymiaru większego niż 1.

5.1.1. Metoda rozkładów warunkowych

Jest to właściwie jedyna metoda ,,w zasadniczy sposób wielowymiarowa”. Opiera się na przedstawieniu gęstości łącznej zmiennych losowych X1,…,Xd jako iloczynu gęstości brzegowej i gęstości warunkowych (wzór łańcuchowy):

f(x1,x2,…,xd)=f(x1)f(x2|x1)f(x3|x1,x2)⋯f(xd|x1,…,xd-1).

Wynika stąd następujący algorytm:

for i:=1 to d do

Gen Xi∼f(⋅|X1,…,Xi-1)

Przykład 5.1 (Wielowymiarowy rozkład normalny)

Ograniczymy się do zmiennych losowych X1,X2 pochodzących z rozkładu dwuwymiarowego N⁢0,0,σ12,σ22,ρ o gęstości

f⁢x1,x2=12⁢π⁢σ1⁢σ2⁢1-ρ2⁢exp⁡-12⁢1-ρ2⁢x12σ12-2⁢ρ⁢x1⁢x2σ1⁢σ2+x22σ22.

Jak wiadomo (można to sprawdzić elementarnym rachunkiem),

f⁢x1=12⁢π⁢σ1⁢exp⁡-x122⁢σ12,

f(x2|x1)=12⁢π⁢σ2⁢1-ρ2exp[-12⁢σ22⁢1-ρ2(x2-ρσ2σ1x1)2].

To znaczy, że N⁢0,σ12 jest rozkładem brzegowym X1 oraz

N⁢ρ⁢σ2σ1⁢x1,σ22⁢1-ρ2

jest rozkładem warunkowym X2 dla X1=x1. Algorytm jst więc następujący.

Gen X1,X2∼N⁢0,1

X1:=σ1⁢X1

X2:=ρ⁢σ2/σ1⁢X1+σ2⁢1-ρ2⁢X2

Rys. 5.1. Próbka z dwuwymiarowego rozkładu normalnego, poziomice gęstości i funkcja regresji x2=E(X2|X1=x1).

Efekt działania powyższego algorytmu widać na Rysunku 5.1. W tym konkretnym przykładzie X1 ma rozkład brzegowy N⁢0,1, zaś X2 ma rozkład warunkowy N⁢X1,0.5. Zauważmy, że Var⁢X2=1.5 i Cov⁢X1,X2=1. Prosta x2=x1 jest wykresem funkcji regresji E(X2|X1=x1) i jest przedstawiona na wykresie. Pokazane są też poziomice gęstości (elipsy). Warto zwrócić uwagę, że funkcja regresji nie pokrywa się ze wspólną osią tych elips. Dlaczego? Jak obliczyć oś?

Uogólnienie na przypadek rozkładu normalnego N⁢μ1,μ2,σ12,σ22,ρ, z niezerowymi średnimi, jest banalne. Uogólnienie na wyższe wymiary też nie jest skomplikowane. Okazuje się jednak, że metoda rozkładów warunkowych dla rozkładów normalnych prowadzi to do algorytmu identycznego jak otrzymany metodą przekształceń.

5.1.2. Metoda przeksztalceń

Podstawą jest następujące twierdzenie o przekształcaniu gęstości.

Twierdzenie 5.1

Załóżmy, że X jest d-wymiarową zmienną losową o wartościach w otwartym zbiorze A⊆Rd. Jeżeli

X∼f⁢x

i h:A→Rd jest dyfeomorfizmem, to

Y=h⁢X∼g⁢y=f⁢h-1⁢y⁢det⁡∂∂⁡y⁢h-1⁢y.

Przykład 5.2 (Wielowymiarowy rozkład normalny)

Rozważmy niezależne zmienne losowe Z1,…,Zd∼N⁢0,1. Wektor Z=Z1,…,Zd ma d-wymiarowy rozkład normalny N⁢0,I o gęstości

f⁢z=2⁢π-d/2⁢exp⁡-12⁢z⊤⁢z.

Jeżeli teraz R jest nieosobliwą macierzą d×d to przekształcenie z↦x=R⁢z jest dyfeomorfizmem z jakobianem det⁡R. Z Twierdzenia 5.1 wynika, że wektor losowy X=R⁢Z ma rozkład normalny o gęstości

fX⁢x=2⁢π-d/2⁢det⁡R-1/2⁢exp⁡-12⁢x⊤⁢Σ-1⁢x,

gdzie Σ=RR⊤. Innymi słowy, X∼N⁢0,Σ. Algorytm generacji jest oczywisty:

Gen Z∼N⁢0,I

X:=R⁢Z

Jeśli dana jest macierz kowariancji Σ wektora X, to przed uruchomieniem algorytmu trzeba znaleźć taką macierz R, żeby Σ=RR⊤. Istnieje wiele takich macierzy, ale najlepiej skorzystać z rozkładu Choleskiego i wybrać macirz trójkątną.

5.2. Kilka ważnych przykładów

Rozmaitych rozkładów wielowymiarowych jest więcej, niż gwiazd na niebie – a wśród nich tyle przykładów ciekawych i ważnych! Musiałem wybrać zaledwie kilka z nich. Rzecz jasna, wybrałem te, które mi się szczególnie podobają.

5.2.1. Rozkłady sferyczne i eliptyczne

Najważniejszy, być może, przykład to wielowymiarowy rozkład normalny, omówiony już w 5.2. Rozpatrzymy teraz rozklady jednostajne na kuli

Bd=x∈Rd:x2≤1

i sferze

Sd-1=x∈Rd:x2=1.

Oczywiście, x oznacza normę euklidesową, x=x12+⋯+xd21/2=x⊤⁢x1/2. Rozkład U⁢Bd ma po prostu stałą gęstość względem d-wymiarowej miary Lebesgue'a na kuli. Rozkład U⁢Sd-1 ma stałą gęstość względem d-1-wymiarowej miary ,,powierzchniowej” na sferze. Oba te rozkłady są niezmiennicze względem obrotów (liniowych przekształceń ortogonalnych) Rd. Takie rozkłady nazywamy sferycznie symetrycznymi lub krócej: sferycznymi. Zauważmy, że zmienną losową o rozkładzie U⁢Sd-1 możemy interpretować jako losowo wybrany kierunek w przestrzeni d-1-wymiarowej. Algorytmy ,,poszukiwań losowych” często wymagają generowania takich losowych kierunków.

Rozkłady jednostajne na kuli i sferze są blisko ze sobą związane.

Jeśli V=V1,…,Vd∼U⁢Bd i R=V to

Y=VR=V1R,…,VdR∼U⁢Sd-1.

Łatwo też zauważyć, że R jest zmienną losową o rozkładzie PR≤r=rd niezależną od Y.
Jeśli Y∼U⁢Sd-1 i R jest niezależną zmienną losową o rozkładzie PR≤r=rd to V=R⁢Y=R⁢Y1,…,R⁢Vd∼U⁢Bd.

Zmienną R łatwo wygenerować metodą odwracania dystrybuanty.

Najprostszy wydaje się algorytm eliminacji:

repeat

Gen V1,…,Vd∼U⁢0,1

until R2=V12+…+Vd2≤1

Na wyjściu otrzymujemy, zgodnie z żądaniem

V=V1,…,Vd∼U⁢Bd.

W istocie, dokładnie ta metoda, dla d=2, była częścią algorytmu biegunowego Marsaglii. Problem w tym, że w wyższych wymiarach efektywnosc eliminacji gwałtownie maleje. Prawdopodobieństwo akceptacji jest równe stosunkowi ,,objętości” kuli Bd do kostki -1,1d. Ze znanego wzoru na objętość kuli d-wymiarowej wynika, że

Bd2d=2⁢πd/2d⁢Γ⁢d/2⋅12d=πd/2d⁢2d-1⁢Γ⁢d/2⟶d→∞0.

Zbieżność do zera jest bardzo szybka. Dla dużego d kula jest znikomą częścią opisanej na niej kostki.

Inna metoda, którą z powodzeniem stosuje się dla d=2 jest związana ze współrzędnymi biegunowymi:

Gen Φ∼U⁢0,2⁢π;

Y1:=cos⁡Φ;⁢Y2:=sin⁡Φ;

Gen U;⁢R:=U;

V1:=Y1⋅R;⁢V2:=Y2⋅R;

Na wyjściu Y1,Y2∼S1 i V1,V2∼B2. Jest to część algorytmu Boxa-Müllera. Uogólnienie na przypadek d>2 nie jest jednak ani proste, ani efektywne. Mechaniczne zastąpienie, współrzędnych biegunowych przez współrzędne sferyczne (dla, powiedzmy d=3) prowadzi do niepoprawnych wyników. Popatrzmy na punkty produkowane przez następujący algorytm:

Gen Φ∼U⁢0,2⁢π;⁢Ψ∼U⁢-π2,π2

Y1:=cos⁡Φ⁢cos⁡Ψ;⁢Y2:=sin⁡Φ⁢cos⁡Ψ;⁢⁢Y3=sin⁡Ψ

$Algorytm ,,współrzędnych sferycznych'' \text@underline{nie generuje} rozkładu jednostajnego na sferze$

Rys. 5.2. Niepoprawne i poprawne generowanie z rozkładu jednostajnego na sferze.

Efekt działania algorytmu ,,współrzędnych sferycznych” jest widoczny na trzech początkowych obrazkach na Rysunku 5.2. Górny lewy rysunek przedstawia punkty widoczne ,,z płaszczyzny równika”, czyli Y1,Y3. Górny prawy – te same punkty widziane ,,znad bieguna”, czyli Y1,Y2. Dolny lewy – górną półsferę widzianą ,,skośnie”, czyli Y2,Y3a, gdzie Y3a=2⁢Y3+2⁢Y1. Dla porównania, na ostatnim rysunku po prawej stronie u dołu – punkty wylosowane rzeczywiście z rozkładu U⁢S2 przy pomocy poprawnego algorytmu podanego niżej.

Gen Z1,…,Zd∼N⁢0,1;

R:=Z12+⋯+Zd21/2;

Y1:=Z1/R,…,Yd:=Zd/R;

Gen U;⁢R:=U1/d;

V1:=Y1⋅R,…,Vd:=Yd⋅R;

Na wyjściu Y∼U⁢Sd-1 i V∼U⁢Bd. Jak widać, algorytm polega na normowaniu punktów wylosowanych ze sferycznie symetrycznego rozkładu normalnego. Jest prosty, efektywny i godny polecenia.

Przykład 5.3 (Wielowymiarowe rozkłady Studenta)

Niech Z=Z1,…,Zd⊤ będzie wektorem losowym o rozkładzie N⁢0,I, zaś R2 – niezależną zmienną losową o rozkładzie χ2⁢n. Wektor

Y1,…,Yd⊤=Z1,…,Zd⊤R2/n

ma, z definicji, Sferyczny rozkład t-Studenta z n stopniami swobody. Gęstość tego rozkładu (z dokładnością do stałej normującej) jest równa

f⁢y=f⁢y1,…,yd∝1+1n⁢∑yi2-n+d/2=1+1n⁢y2-n+d/2.

W przypadku jednowymiarowym, a więc przyjmując d=1, otrzymujemy dobrze znane rozklady t-Studenta z n stopniami swobody o gęstości

f⁢y∝11+y2/nn+1/2

W szczególnym przypadku, biorąc za liczbę stopni swobody n=1, otrzymujemy rozkłady Cauchy'ego. Na przykład, dwuwymiarowy rozkład Cauchy'ego ma taką gęstość:

f⁢y1,y2∝11+y12+y223/2.

Użytecznym uogólnieniem rozkładów sferycznych są rozkłady eliptyczne. Są one określone w następujący sposób. Niech Σ będzie macierzą symetryczną i nieosbliwą. Nazwijmy uogólnionym obrotem przekształcenie liniowe, które zachowuje uogólnioną normę xΣ-1=x⊤⁢Σ-1⁢x1/2. Rozkład jest z definicji eliptycznie konturowany lub krócej eliptyczny, gdy jest niezmienniczy względem uogólnionych obrotów (dla ustalonej macierzy Σ).

5.2.2. Rozklady Dirichleta

Definicja 5.1

Mówimy, że n-wymiarowa zmienna losowa X ma rozkład Dirichleta,

X=X1,…,Xn∼Dir⁢α1,…,αn

jesli X1+⋯+Xn=1 i zmienne X1,…,Xn-1 mają gęstość

f⁢x1,…,xn-1=Γ⁢α1+⋯+αnΓ⁢α1⁢…⁢Γ⁢αn⁢x1α1-1⁢…⁢xn-1αn-1-1⁢1-x1-…⁢xn-1αn-1.

Parametry α1,…,αn mogą być dowolnymi liczbami dodatnimi.

Uwaga 5.1

Rozkłady Dirichleta dla d=2 są to w istocie rozkłady beta:

X1∼Beta⁢α1,α2⁢wtedy i tylko wtedy, gdy ⁢X1,X2∼Dir⁢α1,α2

Wniosek 5.1

Jeśli U1,…,Un-1 są niezależnymi zmiennymi o jednakowym rozkładzie jednostajnym U⁢0,1 i

U1:n<…<Un-1:n-1

oznaczają statystyki pozycyjne, to s⁢p⁢a⁢c⁢j⁢e

Xi=Ui:n-Ui-1:n,⁢Xn=1-Un-1:n

maja rozklad Dir⁢1,…,1.

Twierdzenie 5.2

Jeśli Y1,…,Yn są niezależnymi zmiennymi losowymi o rozkładach gamma,

Yi∼Gamma⁢αi

i S=Y1+⋯+Yn to

X1,…,Xn=Y1S,…,YnS∼Dir⁢α1,…,αn.

Wektor losowy X jest niezależny od S.

Dowód

Obliczymy łączną gęstość zmiennych losowych S,X1,…,Xn-1. Ze wzoru na przekształcenie gęstości wynika, że

⁢⁢fS,X1,…,Xn-1⁢s,x1,…,xn-1=fY1,…,Yn⁢x1⁢s,…,xn⁢s⁢⁢∝x1⁢sα1-1⁢e-x1⁢s⁢…⁢xn⁢sαn-1⁢e-xn⁢s⁢∂⁡y1,…,yn∂⁡s,x1,…⁢xn-1⁢⁢∝x1α1-1⁢…⁢xnαn-1⁢sα1+⋯+αn-1⁢e-s,

ponieważ jakobian przekształcenia odwrotnego jest równy sn-1. Wystarczy teraz zauważyć, że

x1α1-1⁢⋯⁢xnαn-1=Dir,

sα1+⋯+αn-1⁢e-s=Gamma.

∎

Wniosek 5.2

Dla niezależnych zmiennych losowych o jednakowym rozkładzie wykładniczym,

Y1,…,Yn∼Ex⁢1,

jeśli S=Y1+⋯+Yn to

X1,…,Xn=Y1S,…,YnS∼Dir⁢1,…,1

Z 5.1 i 5.2 wynika, że następujące dwa algorytmy:

Gen U1,…,Un-1;

Sort U1,…,Un-1;⁢U0=0;⁢Un=1;

Xi:=Ui-Ui-1

oraz

Gen Y1,…,Yn∼Ex⁢1

S:=Y1+⋯+Yn

Xi:=Yi/S

dają te same wyniki.

Wiele ciekawych własności rozkładów Dirichleta wynika niemal natychmiast z 5.2 (choć nie tak łatwo wyprowadzić je posługując się gęstością 5.2). Mam na myśli przede wszystkim zasadniczą własność ,,grupowania zmiennych”.

Wniosek 5.3

Rozważmy rozbicie zbioru indeksów na sumę rozłącznych podzbiorów:

1,…,n=⋃j=1kIj.

Jeżeli X1,…,Xn∼Dir⁢α1,…,αn i rozważymy ,,zgrupowane zmienne”

Sj=∑i∈IjXi,

to wektor tych zmiennych ma też rozkład Dirichleta,

S1,…,Sk∼Dir⁢β1,…,βk,

gdzie

βj=∑i∈Ijαi.

Co więcej, każdy z wektorów Xi/Sji∈Ij ma rozkład Dirichleta Dir⁢αii∈Ij i wszystkie te wektory są niezależne od S1,…,Sk.

Przykład 5.4

Wniosek 5.3 razem z 5.1 pozwala szybko generować wybrane statystyki pozycyjne. Na przykład łączny rozkład dwóch statystyk pozycyjnych z rozkładu jednostajnego jest wyznaczony przez rozkład trzech ,,zgrupowanych spacji”:

Uk:n-1,⁢Ul:n-1-Uk:n-1, 1-Ul:n-1∼Dir⁢k,l-k,n-l

Przykład 5.5 (Rozkład dwumianowy)

Aby wygenerować zmienną o rozkładzie dwumianowym Bin⁢n,p wystarczy rozpoznać między którymi statystykami pozycyjnymi z rozkładu U⁢0,1 leży liczba p. Nie musimy w tym celu generować wszystkich staystyk pozycyjnych, możemy wybierać ,,najbardziej prawdopodobne”. W połączeniu 5.4 daje to następujący algorytm.

k:=n; θ:=p; X:=0;

repeat

i:=1+k⁢θ;

Gen V∼Beta⁢i,k+1-i;

if θ<V then

begin θ:=θ/V; k:=i-1 end

else

begin X:=X+i; θ=θ-V/1-V; k:=k-i end

until k=0

Metoda generowania zmiennych o rozkładzie Dirichleta opiera sie na następującym fakcie, który jest w istocie szczególnym przypadkiem ,,reguły grupowania” 5.3.

Wniosek 5.4

Jeśli X1,…,Xn∼Dir⁢α1,…,αn i, dla k=1,…,n, określimy Sk=X1+⋯+Xk to zmienne

Y1=S1S2,Y2=S2S3,…,Yn=Sn-1Sn

są niezależne i

Yk∼Beta⁢α1+⋯+αk,αk+1.

Odwrotnie, jeśli zmienne Y1,…,Yn są niezależne i każda z nich ma rozkład beta, to wektor X1,…,Xn ma rozkład Dirichleta.

Oczywiście, jeśli wygenerujemy niezależne zmienne Y1,…,Yn o rozkładach beta, to zmienne X1,…,Xn łatwo ,,odzyskać” przy pomocy wzorów:

Xn=1-Yn-1

Zn-1=1-Yn-2⁢Yn-1

…

Z2=1-Y1⁢Y2⁢⋯⁢Yn-1

Z1=Y1⁢Y2⁢⋯⁢Yn-1

Powyższe równania określają algorytm generowania zmiennych o rozkładzie Dir⁢α1,…,αn.

Treść automatycznie generowana z plików źródłowych LaTeXa za pomocą oprogramowania wykorzystującego LaTeXML.

Symulacje stochastyczne i metody Monte Carlo