Véletlenszámgenerálás elutasításos módszerrel Gamma-eloszlás esetén

 

A segédprogramot Sass Enikõ írta.


A matematikai háttér:

Az elutasításos módszer egy hatékony, általánosan használt eljárás véletlen eltérések generálására, melyek p(x)dx sûrûségfüggvénye (egy érték x és x+dx közötti elõfordulásának valószínûsége) ismert és kiszámítható.

Rajzoljuk meg a grafikonját a p(x) sûrûségfüggvénynek, melyet generálni akarunk, így x bármely tartományában a görbe alatti terület megfelel az adott x elõállításához tartozó valószínûségnek. Ha lenne olyan módszerünk, mellyel a görbe alatti területen (két dimenzióban) egyforma valószínûséggel kiválasztanánk egy véletlen pontot, akkor ezen pont x értéke a kívánt eloszlást mutatná.
Most ugyanerre az ábrára rajzoljunk egy másik, f(x) görbét, mely véges területet határol és mindenütt az eredeti valószínûségi eloszlás fölé esik. (Ez mindig lehetséges, mivel az eredeti görbe, a valószínûség definíciója szerint csak egységnyi területet zár be.) Ezt az f(x)-et összehasonlító függvénynek nevezzük. Tegyük fel, hogy rendelkezünk olyan módszerrel, mellyel egyenlõ valószínûséggel határozhatunk meg véletlen pontot az összehasonlító függvény területe alatt (két dimenzióban). Amikor ez a pont az eredeti valószínûségi eloszlás alatti területen kívül esik, akkor elutasítjuk és egy másik véletlen pontot választunk. Amikor a pont az eredeti valószínûségi eloszlás alatti területre esik, akkor elfogadjuk. Nyilvánvaló, hogy az elfogadott pontok egyenletesek az említett területen, így x értékük a kívánt eloszlással rendelkezik. Látható, hogy a visszautasított pontok hányada csupán a két függvény alatti terület arányán múlik és nem a függvények alakján. Az elfogadott és visszautasított pontok aránya megegyezik a p alatti terület, valamint a p és f közötti terület arányával.
Hogyan válasszunk az f(x) összehasonlító függvény területe alatt (két dimenzióban) egyenletes véletlen számot? Bizonyosodjunk meg róla, hogy az általunk választott összehasonlító függvény határozatlan integrálja analitikusan ismert, valamint analitikusan át is alakítható, hogy megadja x- et "az összehasonlító függvény alatti terület x-tõl balra esõ része" függvényeként. Most válasszunk egy 0 és A közötti egyenletes eltérést, ahol A az f(x) alatti teljes terület, hogy megfelelõ x-et kapjunk. Ezután válasszunk 0 és f(x) között is egyet, mint a kétdimenziós pont y értéke. Belátható, hogy az (x,y) pont egyenletes eloszlású az f(x) összehasonlító függvény alatt.
Természetesen ezt az eljárást átlagosan A-szor kell megismételni, mielõtt a végsõ eltérést megkapnánk.

Az a pozitív egészhez tartozó gamma- eloszlás az a-dik esemény bekövetkeztéig eltelt várakozási idõ, Poisson- folyamat esetén. Például, ha a=1 ez éppen az exponenciális- eloszlást adja, azaz az elsõ esemény várakozási idejét.
A gamma- eloszlás megadja az x és x+dx közé esõ érték elõfordulási valószínûségét.
Az a kis értékei esetén az eltérések elõállításához legjobb ha a darab exponenciális- eloszlású várakozási idõt összeadunk, úgymint az egyenletes eltérések logaritmusait.
Az a nagyobb értékei esetén az eloszlás egy "harangformához" közelít, x = a -beli maximummal és körüli félszélességgel.
A gamma- eloszlás sorbanállási problémák megoldásában lehet hasznos. A ferde eloszlással rendelkezõ valószínûségi változók vizsgálatára használják.

 

A szubrutin használata:

A véletlen számot generáló függvény:
FUNCTION gamdev(ia, idum)

Bemenõ adatok:

Kimenõ adat:

A program felhasználja a már említett ran1 véletlenszám- generáló függvényt, mely egy 0 és1 közötti véletlen eltéréssel tér vissza.
FUNCTION ran1(idum)

Bemenõ adat:

 

A szubrutin meghívása:

Függvényt a fõprogramban értékadó utasítással hívunk meg. A gamma eloszlás esetén a gamdev függvényre:
g=gamdev(ia, idum)
Ez a fõprogramban szereplõ sor azt jelenti, hogy a függvény kimenõ értékét kapja meg a g változó. Ez esetben a függvény neve után az aktuális paraméterlista szerepel.

 

Példa a segédprogram használatára:

PROGRAM teszt
INTEGER ia,n,idum,j
REAL g,gamdev
write(*,*)'Kerem ia erteket! ia:'
read(*,*)ia
write(*,*)'Kerem n erteket! n:'
read(*,*)n

...

g=gamdev(ia,idum) <- itt értékadó utasítással hívjuk meg
write(*,*)g
stop
end

 

próbafuttatás

ia értéke: 2

n értéke: 2

Ha jól gépelted ezt kell kapnod.

kapott eredmény: 3,26... és 1,982... (egymás alatt)

 

 

A segédprogramot innen töltheted le: gamma.for



[Vissza a fõoldalra]