Notice: Undefined index: mode in /home/misc/mst/public_html/common.php on line 63 Notice: Undefined index: mode in /home/misc/mst/public_html/common.php on line 65 Notice: Undefined index: mode in /home/misc/mst/public_html/common.php on line 67 Notice: Undefined index: mode in /home/misc/mst/public_html/common.php on line 69 Notice: Undefined variable: base in /home/misc/mst/public_html/lecture.php on line 36 Symulacje stochastyczne – 4. Generowanie zmiennych losowych II. – MIM UW

Zagadnienia

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 UI:=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:=0R:=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
  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:=1i:=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
  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:=Yreturn 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:=Yreturn 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.

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.