Zagadnienia

4.1 Rozkłady dyskretne
4.2 Schematy kombinatoryczne
- 4.2.1 Pobieranie próbki bez zwracania
- 4.2.2 Permutacje losowe
4.3 Specjalne metody eliminacji
- 4.3.1 Iloraz zmiennych równomiernych
- 4.3.2 Gęstości przedstawione szeregami

4. Generowanie zmiennych losowych II. Specjalne metody

4.1. Rozkłady dyskretne

Jak wylosować zmienną losową $I$ o rozkładzie $\mathbb{P}(I=i)=p_{i}$ mając dane $p_{1},p_{2},\ldots$ ? Metoda odwracanie dystrybuanty w przypadku rozkładów dyskretnych (Przykład 3.6) przyjmuje następującą postać. Obliczamy $s_{1},s_{2},\ldots$ , gdzie $s_{k}=\sum _{{i=1}}^{{k}}p_{i}$ . Losujemy $U\sim{\rm U}(0,1)$ i szukamy przedziału $]s_{{I-1}},s_{I}]$ w którym leży $U$ .

Przykład 4.1 (Algorytm prymitywny)

W zasadzie można to zrobić tak:

Gen $U$ , $I:=1$

while $s_{I}\leq U$ do $I:=I+1$ ;

return $I$

Problemem jest mała efektywność tego algorytmu. Jedno z możliwych ulepszeń polega na bardziej ,,inteligentnym” lokalizowaniu przedziału $]s_{{I-1}},s_{I}]\ni U$ , na przykład metodą bisekcji.

Przykład 4.2 (Algorytm podziału odcinka)

Załóżmy, że mamy rozkład prawdopodobieństwa $(p_{i})$ na przestrzeni skończonej i liczby $p_{i}$ są uporządkowane: $p_{1}>\ldots>p_{m}$ .

Gen $U$ ;

$L:=0$ ; $R:=m$ ;

repeat

$I:=\lfloor({L+R})/{2}\rfloor$ ;

if $U>s_{I}$ then $L:=I$ else $R:=I$ ;

until $L\geq R-1$

return $I$

Przedstawię teraz piękną metodę opartą na innym pomyśle.

Przykład 4.3 (Metoda ,,Alias”)

Przypuśćmy, że mamy dwa ciągi liczb: $q_{1},\ldots,q_{m}$ , gdzie $0<q_{i}<1$ oraz $a(1),\ldots,a(m)$ , gdzie $a(i)\in\{ 1,\ldots,m\}$ (to są owe ,,aliasy”). Rozaptrzmy taki algorytm:

Gen $K\sim{\rm U}\{ 1,\ldots,M\}$ ;

Gen $U$ ;

if $U<q_{K}$ then $I:=K$ else $I:=a(K)$ ;

return $I$

Jest jasne, że

$\mathbb{P}(I=i)=\frac{1}{m}\left(q_{i}+\sum\limits _{{j:a(j)=i}}(1-q_{j})\right).$

Jeśli mamy zadany rozkład prawdopodobieństwa $(p_{1},\ldots,p_{m})$ i chcemy, żeby $\mathbb{P}(I=i)=p_{i}$ , to musimy dobrać odpowiednio $q_{i}$ i $a(i)$ . Można zawsze tak zrobić, i to na wiele różnych sposobów. Opracowanie algorytmu dobierania wektorów $q$ i $a$ do zadanego $p$ pozostawiamy jako ćwiczenie.

4.2. Schematy kombinatoryczne

4.2.1. Pobieranie próbki bez zwracania

Spośród $r$ obiektów chcemy wybrać losowo $n$ tak, aby każdy z $\binom{r}{n}$ podzbiorów miał jednakowe prawdopodobieństwo. Oczywiście, można losować tak, jak ze zwracaniem, tylko odrzucać elementy wylosowane powtórnie.

Przykład 4.4 (Losowanie bez zwracania I)

W tablicy $c(1),\ldots,c(r)$ zaznaczamy, które elementy zostały wybrane.

for $i:=1$ to $r$ do $c(i):=$ false;

$i:=0$

repeat

Gen $K\sim{\rm U}\{ 1,\ldots,r\}$

until $c(K)=$ false;

$c(K):=$ true;

$i:=i+1$ ;

until $i=n$

Pewnym ulepszeniem tej prymitywnej metody jest następujący algorytm.

Przykład 4.5 (Losowanie bez zwracania II)

Tablica $c(1),\ldots,c(r)$ ma takie samo znaczenie jak w poprzednim przykładzie. Będziemy teraz ,,przeglądać” elementy $1,\ldots,r$ po kolei, decydując o zaliczeniu do próbki kolejnego elementu zgodnie z odpowiednim prawdopodobieństwem warunkowym. Niech $i$ oznacza liczbę wybranych, zaś $t-1$ - liczbę przejrzanych poprzednio elementów.

for $i:=1$ to $r$ $c(i):=$ false;

$t:=1$ ; $i:=0$

repeat

Gen $U$ ;

if $U\leq\displaystyle\frac{n-i}{r-(t-1)}$ then

begin

$c(t):=$ true;

$i:=i+1$

end;

$t:=t+1$ ;

until $i=n$

Przykład 4.6 (Losowanie bez zwracania III)

Tablica $s(1),\ldots,s(n)$ będzie teraz zawierała numery wybieranych elementów. Podobnie jak poprzednio, $i$ – oznacza liczbę elementów wybranych, $t-1$ – liczbę przejrzanych.

for $i:=1$ to $n$ do $s(i):=i$ ;

for $t:=n+1$ to $r$ do

begin

Gen $U$ ;

if $U\leq{n}/{t}$ then

begin

Gen $K\sim{\rm U}\{ 1,\ldots,n\}$ ;

$s(K):=t$

end

Uzasadnienie poprawności tego algorytmu jest rekurencyjne. Załóżmy, że przed $t$ -tym losowaniem każda próbka wybrana ze zbioru $\{ 1,\ldots,t-1\}$ ma prawdopodobieństwo

$\binom{t-1}{n}^{{-1}}=\frac{n!}{(t-1)\cdots(t-n)}.$

W kroku $t$ ta próbka ,,przeżywa” czyli pozostaje bez zmiany z prawdopodobieństwem $1-n/t$ (jest to prawdopodobieństwo, że próbka wylosowana ze zbioru $\{ 1,\ldots,t\}$ nie zawiera elementu $t$ . Zatem po kroku $t$ każda próbka nie zawierająca elementu $t$ ma prawdopodobieństwo

$\frac{n!}{(t-1)\cdots(t-n)}\cdot\frac{t-n}{t}=\binom{t}{n}^{{-1}}.$

Z prawdopodobieństwem $n/t$ ,,podmieniamy” jeden z elementów próbki (losowo wybrany) na element $t$ . Sprawdzenie, że każda probka zawierająca element $t$ ma po $t$ -tym losowaniu jednakowe prawdopodobieństwo – pozostawiam jako ćwiczenie.

4.2.2. Permutacje losowe

Przez permutację losową rozumiemy uporządkowanie $n$ elementów wygenerowane zgodnie z rozkładem jednostajnym na przestrzeni wszystkich $n!$ możliwych uporządkowań. Permutację liczb $1,\ldots,n$ zapiszemy w tablicy $\sigma(1),\ldots,\sigma(n)$ . Łatwo sprawdzić poprawność następującego algorytmu.

for $i:=1$ to $n$ do $\sigma(i):=i$ ;

for $i:=1$ to $n-1$ do

begin

Gen $J\sim{\rm U}\{ i,i+1,\ldots,n\}$ ;

$Swap(\pi(i),\pi(J))$

end

Funkcja $Swap$ zamienia miejscami elementy $\pi(i)$ i $\pi(J)$ .

4.3. Specjalne metody eliminacji

4.3.1. Iloraz zmiennych równomiernych

Szczególnie często stosowany jest specjalny przypadek metody eliminacji, znany jako algorytm ,,ilorazu zmiennych równomiernych” (Ratio of Uniforms). Zaczniemy od prostego przykładu, który wyjaśni tę nazwę.

Przykład 4.7 (Rozkład Cauchy'ego)

Metodą eliminacji z prostokąta $[0,1]\times[-1,1]$ otrzymujemy zmienną losową $(U,V)$ o rozkładzie jednostajnym na półkolu $\{(u,v):u\geq 0,u^{2}+v^{2}\leq 1\}$ . Kąt $\Phi$ pomiędzy osią poziomą i punktem $(U,V)$ ma, oczywiście, rozkład ${\rm U}(-\pi,\pi)$ .

repeat

Gen $U\sim{\rm U}(0,1)$ ;

Gen $V\sim{\rm U}(-1,1)$

until $U^{2}+V^{2}<1$

$X:={V}/{U}$

Na wyjściu $X\sim{\rm Cauchy}$ , bo

$\mathbb{P}(X\leq x)=\mathbb{P}(V\leq xU)=\frac{1}{\pi}\Phi+\frac{1}{2}=\frac{1}{\pi}\arctan(x)+\frac{1}{2}.$

Ogólny algorytm metody ,,ilorazu równomiernych” oparty jest na następującym fakcie.

Stwierdzenie 4.1

Załóżmy o funkcji $h:\mathbb{R}\to\mathbb{R}$ , że

$h(x)\geq 0,\,\int h(x)dx<\infty.$

Niech zbiór $C_{h}\subset\mathbb{R}^{2}$ będzie określony następująco:

$C_{n}=\left\{(u,v):0\leq u\leq\sqrt{h\left(\frac{u}{v}\right)}\right\},\,\,|C_{n}|<\infty.$

Miara Lebesgue'a (pole) tego zbioru jest skończone, $|C_{h}|<\infty$ , a zatem można mówić o rozkładzie jednostajnym ${\rm U}(C_{h})$ .

Jeżeli $(U,V)\sim{\rm U}(C_{h})$ i $X={V}/{U}$ , to $X$ ma gęstość proporcjonalną do funkcji $h$ ( $X\sim{h}/{\int h}$ ).

Dowód

,,Pole figury” $C_{h}$ jest równe

$\begin{split}|C_{h}|&=\iint\limits _{{C_{h}}}{\rm d}u{\rm d}v=\iint\limits _{{0\leq u\leq\sqrt{h(x)}}}u{\rm d}u{\rm d}x\\ &=\int\limits _{{-\infty}}^{\infty}\int\limits _{0}^{{\sqrt{h(x)}}}{\rm d}u{\rm d}x=\frac{1}{2}\int\limits _{{-\infty}}^{\infty}h(x){\rm d}x<\infty\end{split}$

Dokonaliśmy tu zamiany zmiennych:

$(u,v)\mapsto\left(u,x=\frac{v}{u}\right).$

Jakobian tego przekształcenia jest równy

$\frac{\partial(u,x)}{\partial(u,v)}=\det\begin{pmatrix}1&0\\ {v}/{u^{2}}&{1}/{u}\\ \end{pmatrix}=\frac{1}{u}$

W podobny sposób, na mocy znanego wzoru na gęstość przekształconych zmiennych losowych obliczamy łączną gęstość $(U,X)$ :

$f_{{U,X}}(u,x)=uf_{{U,V}}(u,ux)\quad\text{na zbiorze }\{ 0\leq u\leq\sqrt{h(x)}\}.$

Stąd dostajemy gęstość brzegową $X$ :

$\int\limits _{{0}}^{{\sqrt{h(x)}}}\frac{u\, du}{|C_{h}|}=\frac{h(x)}{2|C_{h}|}.$

∎

Żeby wylosować $(U,V)\sim{\rm U}(C_{h})$ stosuje się zazwyczaj eliminację z prostokąta. Użyteczne jest następujące oszacowanie boków tego prostokąta.

Stwierdzenie 4.2

Jeśli funkcje $h(x)$ i $x^{2}h(x)$ są ograniczone, wtedy

$C_{h}\subseteq[0,a]\times[b_{-},b_{+}],$

gdzie

$a=\sqrt{\sup\limits _{{x}}h(x)},$

$b_{+}=\sqrt{\sup\limits _{{x\geq 0}}[x^{2}h(x)]},\qquad b_{-}=-\sqrt{\sup\limits _{{x\leq 0}}[x^{2}h(x)]}.$

Dowód

Jeśli $(u,v)\in C_{h}$ to oczywiście $0\leq u\leq\sqrt{h({v}/{u})}\leq\sqrt{\sup _{x}h(x)}$ . Załóżmy dodatkowo, że $v\geq 0$ i przejdżmy do zmiennych $(v,x)$ , gdzie $x=v/u$ . Nierówność $u\leq\sqrt{h({v}/{u})}$ jest równoważna $v^{2}\leq x^{2}h(x)$ . Ponieważ $x\geq 0$ , więc dostajemy $v^{2}\leq b_{+}^{2}$ . Dla $v\leq 0$ mamy analogicznie $v^{2}\leq b_{-}^{2}$ .

∎

Ogólny algorytm RU jest następujący:

repeat

Gen $U_{1},U_{2}$ ;

$U:=aU_{1}$ ; $V:=b_{-}+(b_{+}-b_{-})U_{2}$

until $(U,V)\in C_{h}$ ;

$X:=\displaystyle\frac{V}{U}$

Przykład 4.8 (Rozkład normalny)

Niech

$h(x)=\exp\left(-\frac{1}{2}x^{2}\right).$

Wtedy

$C_{h}=\left\{(u,v):0\leq u\leq\exp\left(-\frac{1}{4}\frac{v^{2}}{u^{2}}\right)\right\}=\left\{\frac{v^{2}}{u^{2}}\leq-4\ln u\right\}.$

Zauważmy, że $a=1$ i $b_{+}=-b_{-}=\sqrt{2e^{{-1}}}$ . Otrzymujemy następujący algorytm:

repeat

Gen $U_{1},U_{2}$

$U:=U_{1}$ ; $V:=\sqrt{2{\rm e}^{{-1}}}(2U_{2}-1)$ ;

$X:=\displaystyle\frac{V}{U}$

until $X^{2}\leq-4\ln U$

4.3.2. Gęstości przedstawione szeregami

Ciekawe, że można skonstruować dokładne algorytmy eliminacji bez konieczności dokładnego obliczania docelowej gęstości $f$ . Podstawowy pomysł jest następujący. Niech $f$ , podobnie jak poprzednio, będzie funkcją proporcjonalną do gęstości (nie jest potrzebna stała normująca) i $f\leq g$ . Załóżmy, że mamy dwa ciągi funkcji, przybliżające $f$ z dołu i z góry:

$\underline{f_{n}}\leq f\leq\overline{f_{n}},\qquad\underline{f_{n}}\to f,\:\overline{f_{n}}\to f\quad(n\to\infty).$

Jeśli umiemy ewaluować funkcje $\underline{f_{n}}$ i $\overline{f_{n}}$ to możemy uruchomić następujący algorytm:

repeat

Gen $Y\sim g$ ;

Gen $U$ ;

$W:=Ug(Y)$ ;

repeat

$n:=n+1$ ;

if $W\leq\underline{f_{n}}(Y)$ then

begin $X:=Y$ ; return $X$ ; stop end

until $W>\overline{f_{n}}(Y)$

until FALSE

Zbieżność przybliżeń dolnych i górnych do funkcji $f$ gwarantuje, że ten algorytm na pewno się kiedyś zatrzyma. Fakt, że na wyjściu $X\sim f/\int f$ jest widoczny.

Poniższy algorytm jest znany jako metoda szeregów zbieżnych. Zakładamy, że

$f(x)=\sum _{{i=1}}^{\infty}a_{n}(x),$

przy czym reszty tego szeregu umiemy oszacować z góry przez znane funkcje,

$\left|\sum _{{i=n+1}}^{\infty}a_{n}(x)\right|\leq r_{{n+1}}(x).$

Pseudo-kod algorytmu jest taki:

repeat

Gen $Y\sim g$ ;

Gen $U$ ;

$W:=Ug(Y)$ ;

$s:=0$ ;

$n:=0$ ;

repeat

$n:=n+1$ ;

$s:=s+a_{n}(Y)$ ;

until $|s-W|>r_{{n+1}}(Y)$ ;

until $W\leq s$ ;

return $X$

Rzecz jasna, w tym algorytmie zmienna $s$ przybiera kolejno wartości równe sumom cząściowym szeregu, $s_{n}(x)=\sum _{{i=1}}^{n}a_{i}(x)$ .

Metoda szeregów naprzemiennych jest dostosowana do sytuacji gdy

$f(x)=g(x)\sum _{{i=0}}^{\infty}(-1)^{n}a_{n}(x)=g(x)\left[1-a_{1}(x)+a_{2}(x)-a_{3}(x)\ldots\right],$

gdzie $1=a_{0}(x)\geq a_{1}(x)\geq a_{2}(x)\geq\cdots\geq 0$ . Wiadomo, że w takiej sytuacji sumy częściowe przybliżają sumę szeregu na przemian z nadmiarem i z niedomiarem. Ten fakt wykorzystuje algorytm szeregów naprzemiennych:

repeat

Gen $Y\sim g$ ;

Gen $U$ ;

$W:=Ug(Y)$ ;

$s:=0$ ;

$n:=0$ ;

repeat

$n:=n+1$ ; { $n$ nieparzyste }

$s:=s+a_{n}(Y)$ ;

if $U\geq s$ then begin $X:=Y$ ; return $X$ ; stop end

$n:=n+1$ ; { $n$ parzyste }

$s:=s-a_{n}(Y)$ ;

until $U>s$ ;

until FALSE

Przykład 4.9

Rozkład Kołmogorowa-Smirnowa ma dystrybuantę

$F(x)=\sum _{{n=-\infty}}^{\infty}(-1)^{n}{\rm e}^{{-2n^{2}x^{2}}},\quad(x\geq 0)$

i gęstość

$f(x)=8\sum _{{n=1}}^{\infty}(-1)^{{n+1}}n^{2}x{\rm e}^{{-2n^{2}x^{2}}},\quad(x\geq 0).$

Możemy zastosować metodę szeregów naprzemiennych, kładąc

$g(x)=4x{\rm e}^{{-2x^{2}}}\quad(x\geq 0)$

oraz

$a_{n}(x)=n^{2}x{\rm e}^{{-2(n^{2}-1)x^{2}}}.$

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

Symulacje stochastyczne i metody Monte Carlo

Zagadnienia

4. Generowanie zmiennych losowych II. Specjalne metody

4.1. Rozkłady dyskretne

Przykład 4.1 (Algorytm prymitywny)

Przykład 4.2 (Algorytm podziału odcinka)

Przykład 4.3 (Metoda ,,Alias”)

4.2. Schematy kombinatoryczne

4.2.1. Pobieranie próbki bez zwracania

Przykład 4.4 (Losowanie bez zwracania I)

Przykład 4.5 (Losowanie bez zwracania II)

Przykład 4.6 (Losowanie bez zwracania III)

4.2.2. Permutacje losowe

4.3. Specjalne metody eliminacji

4.3.1. Iloraz zmiennych równomiernych

Przykład 4.7 (Rozkład Cauchy'ego)

Stwierdzenie 4.1

Dowód

Stwierdzenie 4.2

Dowód

Przykład 4.8 (Rozkład normalny)

4.3.2. Gęstości przedstawione szeregami

Przykład 4.9