Zagadnienia

4. Generowanie zmiennych losowych II. Specjalne metody

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 UU0,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 UI:=1
  while sIU 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:=0R:=m;
repeat
  I:=L+R/2;
  if U>sI then L:=I else R:=I;
until LR-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 a1,,am, gdzie ai1,,m (to są owe ,,aliasy”). Rozaptrzmy taki algorytm:

Gen KU1,,M;
Gen U;
if U<qK then I:=K else I:=aK;
return I

Jest jasne, że

PI=i=1mqi+j:aj=i1-qj.

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

for i:=1 to r do ci:= false;
i:=0
repeat
  repeat
    Gen KU1,,r
  until cK= false;
  cK:= 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 c1,,cr 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 ci:= false;
t:=1i:=0
repeat
  Gen U;
  if Un-ir-t-1 then
    begin
    ct:= true;
    i:=i+1
  end;
  t:=t+1;
until i=n
Przykład 4.6 (Losowanie bez zwracania III)

Tablica s1,,sn 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 si:=i;
for t:=n+1 to r do
  begin
  Gen U;
  if Un/t then
    begin
    Gen KU1,,n;
    sK:=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,,t-1 ma prawdopodobieństwo

t-1n-1=n!t-1t-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-1t-nt-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 JUi,i+1,,n;
  Swapπi,πJ
  end
Funkcja Swap 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:u0,u2+v21. Kąt Φ pomiędzy osią poziomą i punktem U,V ma, oczywiście, rozkład U-π,π.

repeat
  Gen UU0,1;
  Gen VU-1,1
until U2+V2<1
X:=V/U

Na wyjściu XCauchy, bo

PXx=PVxU=1πΦ+12=1πarctanx+12.

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

Stwierdzenie 4.1

Załóżmy o funkcji h:RR, że

hx0,hxdx<.

Niech zbiór ChR2 będzie określony następująco:

Cn=u,v:0uhuv,Cn<.

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

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

Dowód

,,Pole figury” Ch jest równe

Ch=Chdudv=0uhxududx=-0hxdudx=12-hxdx<

Dokonaliśmy tu zamiany zmiennych:

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

Jakobian tego przekształcenia jest równy

u,xu,v=det10v/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,Xu,x=ufU,Vu,uxna zbiorze 0uhx.

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

0hxuduCh=hx2Ch.

Żeby wylosować U,VUCh 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 hx i x2hx są ograniczone, wtedy

Ch0,a×b-,b+,

gdzie

a=supxhx,
b+=supx0x2hx,b-=-supx0x2hx.
Dowód

Jeśli u,vCh to oczywiście 0uhv/usupxhx. Załóżmy dodatkowo, że v0 i przejdżmy do zmiennych v,x, gdzie x=v/u. Nierówność uhv/u jest równoważna v2x2hx. Ponieważ x0, więc dostajemy v2b+2. Dla v0 mamy analogicznie v2b-2.

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

repeat
  Gen U1,U2;
  U:=aU1V:=b-+b+-b-U2
until U,VCh;
X:=VU

Przykład 4.8 (Rozkład normalny)

Niech

hx=exp-12x2.

Wtedy

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

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

repeat
  Gen U1,U2
  U:=U1;  V:=2e-12U2-1;
  X:=VU
until X2-4lnU

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

fn¯ffn¯,fn¯f,fn¯fn.

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

repeat
Gen Yg;
Gen U;
W:=UgY;
  repeat
  n:=n+1;
  if Wfn¯Y then
     begin X:=Yreturn 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 Xf/f jest widoczny.

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

fx=i=1anx,

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

i=n+1anxrn+1x.

Pseudo-kod algorytmu jest taki:

repeat
Gen Yg;
Gen U;
W:=UgY;
s:=0;
n:=0;
  repeat
  n:=n+1;
  s:=s+anY;
  until s-W>rn+1Y;
until Ws;
return X
Rzecz jasna, w tym algorytmie zmienna s przybiera kolejno wartości równe sumom cząściowym szeregu, snx=i=1naix.

Metoda szeregów naprzemiennych jest dostosowana do sytuacji gdy

fx=gxi=0-1nanx=gx1-a1x+a2x-a3x,

gdzie 1=a0xa1xa2x0. 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 Yg;
Gen U;
W:=UgY;
s:=0;
n:=0;
  repeat
  n:=n+1;     { n nieparzyste }
  s:=s+anY;
  if Us then begin X:=Yreturn X; stop end
  n:=n+1;  { n parzyste }
  s:=s-anY;
  until U>s;
until FALSE

Przykład 4.9

Rozkład Kołmogorowa-Smirnowa ma dystrybuantę

Fx=n=--1ne-2n2x2,x0

i gęstość

fx=8n=1-1n+1n2xe-2n2x2,x0.

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

gx=4xe-2x2x0

oraz

anx=n2xe-2n2-1x2.

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.