Symulacje stochastyczne – 4. Generowanie zmiennych losowych II.

4.1. Rozkłady dyskretne

Jak wylosować zmienną losową I o rozkładzie PI=i=pi mając dane p1,p2,…? Metoda odwracanie dystrybuanty w przypadku rozkładów dyskretnych (Przykład 3.6) przyjmuje następującą postać. Obliczamy s1,s2,…, gdzie sk=∑i=1kpi. Losujemy U∼U⁢0,1 i szukamy przedziału ]sI-1,sI] w którym leży U.

Przykład 4.1 (Algorytm prymitywny)

W zasadzie można to zrobić tak:

Gen U, I:=1

while sI≤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 ]sI-1,sI]∋U, na przykład metodą bisekcji.

Przykład 4.2 (Algorytm podziału odcinka)

Załóżmy, że mamy rozkład prawdopodobieństwa pi na przestrzeni skończonej i liczby pi są uporządkowane: p1>…>pm.

Gen U;

L:=0; R:=m;

repeat

I:=L+R/2;

if U>sI then L:=I else R:=I;

until L≥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: q1,…,qm, gdzie 0<qi<1 oraz a⁢1,…,a⁢m, gdzie a⁢i∈1,…,m (to są owe ,,aliasy”). Rozaptrzmy taki algorytm:

Gen K∼U⁢1,…,M;

Gen U;

if U<qK then I:=K else I:=a⁢K;

return I

Jest jasne, że

PI=i=1m⁢qi+∑j:a⁢j=i1-qj.

Jeśli mamy zadany rozkład prawdopodobieństwa p1,…,pm i chcemy, żeby PI=i=pi, to musimy dobrać odpowiednio qi 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 rn 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,…,c⁢r zaznaczamy, które elementy zostały wybrane.

for i:=1 to r do c⁢i:= false;

i:=0

repeat

Gen K∼U⁢1,…,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,…,c⁢r ma takie samo znaczenie jak w poprzednim przykładzie. Będziemy teraz ,,przeglądać” elementy 1,…,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≤n-ir-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,…,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≤n/t then

begin

Gen K∼U⁢1,…,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,…,t-1 ma prawdopodobieństwo

t-1n-1=n!t-1⁢⋯⁢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,…,t nie zawiera elementu t. Zatem po kroku t każda próbka nie zawierająca elementu t ma prawdopodobieństwo

n!t-1⁢⋯⁢t-n⋅t-nt=tn-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,…,n zapiszemy w tablicy σ⁢1,…,σ⁢n. Łatwo sprawdzić poprawność następującego algorytmu.

for i:=1 to n do σ⁢i:=i;

for i:=1 to n-1 do

begin

Gen J∼U⁢i,i+1,…,n;

S⁢w⁢a⁢p⁢π⁢i,π⁢J

end

Funkcja S⁢w⁢a⁢p zamienia miejscami elementy π⁢i i π⁢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×-1,1 otrzymujemy zmienną losową U,V o rozkładzie jednostajnym na półkolu u,v:u≥0,u2+v2≤1. Kąt Φ pomiędzy osią poziomą i punktem U,V ma, oczywiście, rozkład U⁢-π,π.

repeat

Gen U∼U⁢0,1;

Gen V∼U⁢-1,1

until U2+V2<1

X:=V/U

Na wyjściu X∼Cauchy, bo

PX≤x=PV≤x⁢U=1π⁢Φ+12=1π⁢arctan⁡x+12.

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

Stwierdzenie 4.1

Załóżmy o funkcji h:R→R, że

h⁢x≥0,⁢∫h⁢x⁢d⁢x<∞.

Niech zbiór Ch⊂R2 będzie określony następująco:

Cn=u,v:0≤u≤h⁢uv,⁢⁢Cn<∞.

Miara Lebesgue'a (pole) tego zbioru jest skończone, Ch<∞, a zatem można mówić o rozkładzie jednostajnym U⁢Ch.

Jeżeli U,V∼U⁢Ch i X=V/U, to X ma gęstość proporcjonalną do funkcji h (X∼h/∫h).

Dowód

,,Pole figury” Ch jest równe

⁢Ch⁢=∬Chd⁢u⁢d⁢v=∬0≤u≤h⁢xu⁢d⁢u⁢d⁢x⁢⁢=∫-∞∞∫0h⁢xd⁢u⁢d⁢x=12⁢∫-∞∞h⁢x⁢d⁢x<∞⁢

Dokonaliśmy tu zamiany zmiennych:

(u,v)↦(u,x=vu).

Jakobian tego przekształcenia jest równy

∂⁡u,x∂⁡u,v=det⁡10v/u21/u=1u

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

fU,X⁢u,x=u⁢fU,V⁢u,u⁢x⁢⁢na zbiorze 0≤u≤h⁢x.

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

∫0h⁢xu⁢⁢d⁢uCh=h⁢x2⁢Ch.

∎

Żeby wylosować U,V∼U⁢Ch 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 x2⁢h⁢x są ograniczone, wtedy

Ch⊆0,a×b-,b+,

gdzie

a=supx⁡h⁢x,

b+=supx≥0⁡x2⁢h⁢x,⁢b-=-supx≤0⁡x2⁢h⁢x.

Dowód

Jeśli u,v∈Ch to oczywiście 0≤u≤h⁢v/u≤supx⁡h⁢x. Załóżmy dodatkowo, że v≥0 i przejdżmy do zmiennych v,x, gdzie x=v/u. Nierówność u≤h⁢v/u jest równoważna v2≤x2⁢h⁢x. Ponieważ x≥0, więc dostajemy v2≤b+2. Dla v≤0 mamy analogicznie v2≤b-2.

∎

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

repeat

Gen U1,U2;

U:=a⁢U1; V:=b-+b+-b-⁢U2

until U,V∈Ch;

X:=VU

Przykład 4.8 (Rozkład normalny)

Niech

h⁢x=exp⁡-12⁢x2.

Wtedy

Ch={(u,v):0≤u≤exp(-14v2u2)}={v2u2≤-4lnu}.

Zauważmy, że a=1 i b+=-b-=2⁢e-1. Otrzymujemy następujący algorytm:

repeat

Gen U1,U2

U:=U1; V:=2⁢e-1⁢2⁢U2-1;

X:=VU

until X2≤-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≤g. Załóżmy, że mamy dwa ciągi funkcji, przybliżające f z dołu i z góry:

fn¯≤f≤fn¯,⁢fn¯→f,⁢fn¯→f⁢n→∞.

Jeśli umiemy ewaluować funkcje fn¯ i fn¯ to możemy uruchomić następujący algorytm:

repeat

Gen Y∼g;

Gen U;

W:=U⁢g⁢Y;

repeat

n:=n+1;

if W≤fn¯⁢Y then

begin X:=Y; return X; stop end

until W>fn¯⁢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∼f/∫f jest widoczny.

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

f⁢x=∑i=1∞an⁢x,

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

∑i=n+1∞an⁢x≤rn+1⁢x.

Pseudo-kod algorytmu jest taki:

repeat

Gen Y∼g;

Gen U;

W:=U⁢g⁢Y;

s:=0;

n:=0;

repeat

n:=n+1;

s:=s+an⁢Y;

until s-W>rn+1⁢Y;

until W≤s;

return X

Rzecz jasna, w tym algorytmie zmienna s przybiera kolejno wartości równe sumom cząściowym szeregu, sn⁢x=∑i=1nai⁢x.

Metoda szeregów naprzemiennych jest dostosowana do sytuacji gdy

f⁢x=g⁢x⁢∑i=0∞-1n⁢an⁢x=g⁢x⁢1-a1⁢x+a2⁢x-a3⁢x⁢…,

gdzie 1=a0⁢x≥a1⁢x≥a2⁢x≥⋯≥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∼g;

Gen U;

W:=U⁢g⁢Y;

s:=0;

n:=0;

repeat

n:=n+1; { n nieparzyste }

s:=s+an⁢Y;

if U≥s then begin X:=Y; return X; stop end

n:=n+1; { n parzyste }

s:=s-an⁢Y;

until U>s;

until FALSE

Przykład 4.9

Rozkład Kołmogorowa-Smirnowa ma dystrybuantę

F⁢x=∑n=-∞∞-1n⁢e-2⁢n2⁢x2,x≥0

i gęstość

f⁢x=8⁢∑n=1∞-1n+1⁢n2⁢x⁢e-2⁢n2⁢x2,x≥0.

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

g⁢x=4⁢x⁢e-2⁢x2⁢x≥0

oraz

an⁢x=n2⁢x⁢e-2⁢n2-1⁢x2.

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