Modellezés és szimuláció

2. A szimuláció matematikai alapjai

A számítógépes szimulációt matematikai szempontból két oldalról kívánjuk megközelíteni. Egyrészt a determinisztikus gép „véletlen generálási” lehetőségeit, módjait vizsgáljuk, amivel a számítógép ezen önmagában véve paradoxnak tűnő „képessége” körül oszlatjuk el a ködöt. Másrészt megtanuljuk a számítógép adta véletlenszámok sorozatát természetutánzásra használni, miközben érthetővé válik, hogy milyen alapon várhatunk el a szimulációs modellünktől, egy véletlen irányította modelltől objektív, a modellezettről „igazakat” állító eredményeket.

2.1. véletlenszámok generálása

A számítógép produkálta véletlenről – úgy gondoljuk – nem eshet szó addig, amíg a véletlen „természetét”, a véletlentől elvárt tulajdonságokat nem tisztáztuk.

A. A véletlenszerű kísérletektől elvárjuk, hogy minden lehetséges kimenetele előbb-utóbb bekövetkezzék.

B. Az eseményre vonatkozó kísérleteinket jellemezze az is, hogy addigi bekövetkezésekből (vagy be nem következésekből) ne következtethessünk a továbbiak kimeneteleire. (Azaz legyenek egymástól függetlenek.)

C. E két kvalitatív elvárásunkat pontosíthatjuk az alábbi kvantitatív állításunkkal: a szimulációnk akkor valósághű, ha a szimulált kísérletek kimeneteleinek relatív gyakorisága kifejezi a hozzátartozó valóságos események valószínűségeit. (Azaz a nem nulla valószínűségű esemény relatív gyakorisága nem nulla, pontosabban akkora, mint maga a valószínűség. Másrészt két vagy több kísérletünk kimeneteleinek együttes relatív gyakorisága egyezzék meg ezek együttes bekövetkezésének valószínűségével.)

(Számítógépről lévén szó azon semmi csodálkozni való nincs, hogy egy számalgoritmust keresünk, amely véletlen események gyanánt véletlenszámot eredményez.)

Kérdés tehát, hogy lehet-e olyan determinisztikus algoritmust konstruálni, amely a fentieket teljesíti? Az igenlő választ történetileg először Neumann János adta meg: a négyzetközép-nek elnevezett módszerével. A lényege a következő:

Induljunk ki egy 2*k jegyű számból! Emeljük négyzetre, majd az így keletkezett 4*k jegyű szám középső 2*k jegyét tekintsük az első véletlenszámnak (itt a véletlen esemény maga a szám kijötte)! Ebből a számból állítjuk elő a fenti módon a másodikat, és így tovább ...

A négyzetközép módszer egy variánsát – a szorzatközép módszert – az alábbi algoritmus valósítja meg! Állítsunk elő kétjegyű álvéletlenszámokat, a kiinduló értéket olvassa be a program! (A négyzetre emelés művelet szerepét töltse be az A*előző szám+B formula!)

Ha R0-nak például a 32 kezdőértéket adjuk, az alábbi sorozatot kapjuk:

40, 49, 59, 70, 82, 95, 9, 15, 21, 28, 36, 44, 53, 63, 74, 86, 99, 14, 20, 27, 35, 43, 52, 62, 73, 85, 98, 13, 19, 26, 33, 41, 50, 60, 71, 83, 96, 10, 16, 22, 29, 37, 46, 55, 65, 76, 88, 2, 7, 13, ...

Célszerű több kezdőérték mellett (sőt más A és B esetére) is álvéletlenszám-sorozatokat generálni! Vegyük észre a módszer jellegzetességeit, hibáit!

1. Gyakoriak az olyan kezdőértékek (R0,A,B), amelyekre nagyon rosszul viselkedik: elfajul a sorozat. (Például állandóan ugyanazt az egy értéket szolgáltatja.)

2. Észrevehetően periodikus, ami a determinisztikusságból fakad. (Elvileg sem tudna 100-nál hosszabb periódusú sorozatokat produkálni, hiszen kétjegyű szám csak 100 van, ugyanabból a számból generálva másodszor is csak ugyanúgy folytathatja, azaz az ismétlődés garantált.)

3. A kapott álvéletlenszám-sorozatokat az alábbi szakaszokra lehet bontani:

[aperiodikus szakasz] [periodikus szakasz] [periodikus szakasz] …

Ebből következik, hogy nem célszerű az aperiodikus szakasz hosszánál több véletlenszámot generálni!

4. Aprólékosabban elemezve a kapott sorozatot újabb kellemetlen sajátosságokat fedezhetünk föl. Például hosszabb, növekvő (méghozzá sajnálatosan elég szisztematikusan növekvő!) szakaszok ismétlődését tapasztalhatjuk.

A módszer mindezen, túlontúl „átlátszó” hibái az eljárás lényegéből fakadnak, mégis megemlítettük történeti érdekessége és egyszerű megvalósíthatósága miatt.

Lineáris kongruencia módszer

• v₀:=tetszőleges egész szám (vehetjük a pl. számítógép órájából, hogy ne legyen megismételhető)
• v_n+1:=(a*v_n+c) mod m

Állítás: Ha m=2^k, a=4*x+1, (c,m)=1 (és m prímosztói a–1-nek is prímosztói) , akkor m lesz a periódushossz.

Megjósolhatóság kérdése:

• A v₀ nem ismert.
• Az a,c nem ismert.
• Az m nem ismert? x_i:=v_i/m valós szám! Ha 0≤v_i<m, akkor 0≤x_i<1!

A gyakorlatban használatos véletlenszám-generátorok többsége a lineáris kongruencia módszer speciális esete:

Adott k+1 db egész szám (a₀,a₁,...,a_k), ezek a kongruencia együtthatói, és k db (R₀,R₁,...,R_k-1) kezdőérték. Az R_n+1-t az előzőkből az R_n+1≡a₀*R_n+a₁*R_n-1+ ... +a_k-1*R_n-k+1+a_k (mod M) formula alapján számoljuk.

Ennek kiszámolásáról csak annyit, hogy a következő véletlenszám a ≡-jel jobboldalán levő kifejezés értékének M-mel való osztási maradéka. Az M modulus és az együtthatók alkalmas választásával megközelíthető a periódus elvi korlátja (M^k).

Az egyes konkrét véletlenszám-generátorokhoz kiszámolható a periódus. Hogy némi fogalmunk legyen ennek nagyságrendjéről, közöljük, hogy például egy R_n+1≡(10^a+1)*R_n+R_n-1 (mod 10^b) formulával definiált véletlenszám-generátor periódushosszát az 1.5*10^b képlettel határozhatjuk meg.

A véletlenszám-előállító algoritmusok további példájaként egy érdekes módszert említünk. A kongruencia módszerből M=2 és az első együttható=1 esetén kapjuk. A továbbiakban az együtthatók is, a véletlenszámok is a 0 vagy az 1 értékek valamelyikét vehetik föl. Úgy értelmezhetjük a módszert, mintha a számokat „bitenként” véletlenül generálná. Az érdekessége az, hogy meglepően jó a periodicitás szempontjából: 2^k-1. (A k a következő bit előállításában résztvevő bitek száma.) Ha m bináris jegy alkotja a véletlenszámunkat, akkor a következőt az előző m-1 db, valamint az újonnan generált véletlen bit alkotja. Ily módon éppen 2^k ilyen szám létezik, tehát csak eggyel több, mint a periódus hossza. (Melyik szám lehet az, amit nem generál?) Továbbiakban jelöljük P-vel a periódus hosszát!

Könnyen beláthatók a módszer és a P ismeretében az alábbi, a használhatóságot jelentő állítások.

1. Ha egy kezdő k bit-variációra P=2^k-1 periódushosszú sorozatot ad, akkor tetszőlegesből indítva (kivéve a csupa 0-ból állót) a periódus hossza továbbra is ugyanazon P marad. (Hiszen a k db 0-t tartalmazó szakaszt kivéve az összes bit-k-ast generálja.)

2. Az aperiodikus szakasz hossza is P.

3. Az m (<k+1) bites bit-variációk száma 2^k-m, a csak 0-t tartalmazó m-esek száma pedig eggyel kevesebb. Ennek következménye, hogy „majdnem” egyenletes eloszlást követ.

Persze ez utóbbi módszer ritka példája annak, hogy a véletlenszám generátornak elvileg határozhatjuk meg az egyes számvariációk „eloszlását”. S ezzel bizonyíthatjuk jóságát vagy éppen használhatatlanságát. A gyakorlatban nem ilyen kedvező az ellenőrizhetőség helyzete. Általában mintasorozatokat generálva készítenek statisztikákat, s mint egy valódi véletlenfolyam tulajdonságaira ily módon következtetnek.

Egy kis kitérő: Vajon melyek a „jóság” legfontosabb kritériumai?

Vegyünk egy v_i (i=1,...) véletlenszám sorozatot!

• a sorozat 1-egyenletes – a v_i véletlenszámok a [0,M) intervallum bármely [a,b) részintervallumába esés valószínűsége csak az intervallum hosszától függ;
• a sorozat 2-egyenletes – a (v_i,v_i+1) véletlenszám párok a ([0,M), [0,M)) négyzet bármely ([a,b),[c,d)) résztéglalapjába esés valószínűsége csak a téglalap területétől függ;
• a sorozat K-egyenletes – ...;
• a sorozat ∞-egyenletes – minden K természetes számra K-egyenletes.

Mivel minden K-ra nem tudjuk megnézni, más lehetőségeket kell választanunk:

• k-hosszúságú sorozatainak (variációinak) gyakorisága (k=1,2, ...),
• adott nagyságrendű számjegyéből képzett k-hosszúságú sorozatainak (variációinak) gyakorisága (k=1,2,...),
• egy szám (vagy egy adott helyiértékű számjegy) két előfordulása közötti „hézagok” vizsgálata,
• statisztikus próbák (pl. illeszkedés- vagy függetlenségvizsgálatra a χ-négyzet próba),
• adott szám- (vagy számjegy-, esetleg bit-) minták gyakoriságai,
• kombinációk gyakoriságai (Póker-teszt),
• az egyenletesség-ellenőrzés egy látványos (szubjektív) módszere a következő: a véletlenszámok sorozatát bontsuk számpárok sorozatára. E számpárokat tekintsük egy síkbeli pont koordinátáinak! Leképezve a képernyőre ezeket a pontokat, gyújtsuk ki a képernyőn! Gondoljuk meg, a véletlenszám-generátor egyenletességét feltételezve: mi mondható a képernyő „kitöltődéséről”, és ennek üteméről!

Az utóbbi módszer elvi alapját Weyl-től származó alábbi tétel alkotja:

Annak, hogy az n-dimenziós kockában véletlenszerűen elhelyezett pontsorozat egyenletes eloszlású legyen, szükséges és elégséges feltétele, hogy bármely T tartomány esetén

ahol Q(T) a teljes kocka térfogatát, H(T,n) a T-be esés relatív gyakoriságát, V(T) pedig a T térfogatát jelöli.

A K-számjegy-kombinációk gyakoriságát számláló algoritmus az alábbi! Gyak(i) -ben számláljuk az i. kombináció gyakoriságát, a Szám() vektorban tároljuk a következő K véletlen számjegyet.

A program kifejtésébe most nem bonyolódunk. Az Olvasónak maradjon gyakorlásként a „specifikált” eljárások részletezése!

Kombinált módszerek

A véletlenszám-generátorok egyfajta javítási lehetőségével ismerkedünk meg a téma lezárásaként:

Adott két (esetleg több) P1, P2 (...) periódushosszú generátor. A két (több) generátor egy (valamilyen „topológiájú”) integrált rendszerben állítja elő a véletlenszámokat. Ennek topológiája szerint a következő felosztásban csoportosíthatjuk a véletlenszám-generátorokat:

• párhuzamosított: Az f egy „egyszerű” függvénye a két véletlenszámnak; például hol az egyiket, hol a másikat adja tovább (azazmintegy összefésüli), esetleg a kettőt eggyé szuperponálja.

Két generátor párhuzamos kapcsolása.

• visszacsatolásos: A következő véletlenszámot a két véletlenszám-generátor „valódi kooperációjával” állítjuk elő úgy, hogy mindegyik a másik előzőjét veszi figyelembe a soron következő véletlenszám generálásához. Az f hasonló, mint az előző esetben. Ennek speciális esete az a módszer, amelyben az R2-nek a „zavar” szerepe jut. Az R1 által generált fő véletlenszám-sorozat minden m. elemét az R2 adja, s így „zavarja” az R1 működését. (Innen kapta a zavarásmódszer elnevezést is.)

Két generátor visszacsatolásos kapcsolása.

• sorosan kapcsolt: Mindkét generátor a másik előző számából állítja elő a következő véletlenszámot.

Két generátor soros kapcsolása.

Sejthető, hogy ezen kombinált módszerrel kapott véletlen-sorozat általában kedvezőbb lesz periodicitási (periódushossz és stabilitás – kezdőértéktől való kisebb függés) és statisztikus szempontokból, de nem jár a megvalósító program lényeges bonyolításával. (Megjegyezzük, hogy a módszerrel óvatosan kell bánni, néha váratlanul rossz sorozatot eredményez.)

A szorzatközép módszerre megírt programot felhasználva adjuk meg két különböző A, B paraméterrel működő generátor „összefésülésének” algoritmusát!

A véletlenszámokkal szemben támasztott „heurisztikus” követelményeink közül a másodikra (B) térjünk vissza még egy, talán meglepőnek tűnő gondolattal! Sokszor még az sem lényeges, hogy legalább „látszólag” megjósolhatatlanul kövessék egymást a véletlenszámok. Olyan esetekre kell gondolnunk, amelyeknél a szimulációs folyamat „egyszerű” struktúrájú. Ilyenkor bármely olyan számsorozat jó, amelyben minden szükséges variáció az elméletileg elvárt gyakorisággal fordul elő. Az előfordulásuk sorrendje lényegtelen! Akár növekvően rendezve is adhatók.

1. példa: ha 1 és 9 közötti véletlenszámokra van szükségünk (és a szimulált rendszer a fenti tulajdonságú), akkor a következő sorozat is megfelelő:

A folytatás – úgy gondoljuk – mindenki számára világos. Próbálja a fenti sorozatot program segítségével generálni! A sorozat felhasználásáról csak annyit, hogy abbahagyni „10^k-1 felbontásánál” (valamely k-ra) szabad. (Ugyanis ekkor garantálható legalább a számjegyek azonos gyakorisága!).

2. példa: ha egy k dimenziós valószínűségi vektorváltozót használunk föl a szimuláció során, akkor k-ad rendig bezárólag kell egyenletesen szerepeltetni a variációkat, így:

Vissza a tartalomjegyzékhez

2.2. A szimuláció elvi alapjai

A számítógépes szimuláció elvi „jóságának” indoklását a sztochasztikus folyamatokkal való kézenfekvő kapcsolata segítségével építjük föl.

A kapcsolat elvi alapjait négy aspektusból fogjuk áttekinteni.

A. A szimuláció és a szimulálandó elvi azonosságát, mint folyamatok azonosságát vizsgáljuk.

B. Látni fogjuk, hogy a szimulációs folyamatot a szimuláció közben „állóképek” sorozatából rakhatjuk össze. Csak ezen pillanatfelvételek jóságát kell biztosítanunk. Ehhez annak vizsgálata szükséges, hogy a folyamat adott, rögzített időpontbeli állapoteloszlását miként tudjuk valósághűen visszatükrözni.

C. Az elvi alapok taglalását a szimuláció „módszertani” alapjául szolgáló Glivenko-tétellel folytatjuk.

D. Végül egy a szimuláció praktikus vetülete felé átvezető problémát vetünk föl: meddig kell folytatnunk a kísérletet, hogy „megnyugtatónak” érezhessük az eredményeinket.

Az elméleti alapok ismertetése után rátérünk a számítógépes szimuláció gyakorlati kérdéseire (a 3.3. pontban). Itt a számítógép véges pontosságából adódó hiba vizsgálatával kezdjük a gyakorlati problémák sorát, majd számos programozástechnikai részlettel, fogással folytatjuk.

A. A szimulációs folyamat, mint „filmkockák sorozata”

Írja le a szimulálandó folyamatot az X(t), a szimulációs folyamatot az Y(t) valószínűségi változó. A szimulációt akkor tekinthetjük jónak, ha a szimulálandó és a szimulációs folyamat (bármely t időpillanatában) egy tetszőlegesen választott állapotában azonos valószínűséggel tartózkodik.

Azaz, P(X(t)=s)=P(Y(t)=s), feltéve, hogy az „indulás” körülményei azonosak voltak: P(X(0)=r)=P(Y(0)=r), ahol r és s tetszőleges folyamatbeli állapot. Igaz a következő állítás:

A P(X(t)=s)=P(Y(t)=s) teljesülésének szükséges és elégséges feltétele, hogy a P(X(t)=s|∀u<t: X(u)∈Γ) = P(Y(t)=s|∀u<t: Y(u)∈Γ) ∀s állapotra. Itt a Γ a folyamat pályáját jelöli.

Az állítás jelentése a következő: sztochasztikus értelemben a szimulációs folyamat akkor és csak akkor ekvivalens a szimulálandóval, ha a szimuláció során az állapotváltozás valószínűségi szempontból azonos a szimulálandóéval, figyelembe véve a folyamat eddig befutott pályáját. Tehát a folyamat teljes ismerete helyett (hiszen általában ennek megismerése a cél) elegendő a kezdeti állapoteloszlás valamint az állapotváltozás módjának ismerete, ami a legrosszabb esetben is a teljes eddigi pálya figyelembe vételét jelenti. A konkrét esetekben látni fogjuk, hogy nagyon sokszor ennél jóval kevesebb információval is meghatározható a továbblépés eloszlása.

B. A folyamat „filmkockáiról”

Figyelmünket most a folyamat egy adott, rögzített időpontbeli „pillanatképére” fordítjuk. Induljunk ki abból a feltevésből, hogy rendelkezünk a szimulációhoz egy F eloszlást követő X valószínűségi változóval. Kérdés, hogyan állíthatunk elő ebből más eloszlású valószínűségi változókat? Erre egy direkt formulával válaszolunk, amely alkalmazására később számos példát látunk. Ha a G eloszlású Y valószínűségi változót kell szimulálnunk, akkor az X-et transzformálnunk kell, mégpedig a G^-1(F(X)) formula segítségével. (Itt a G^-1 a G inverzét jelöli!)

E transzformáció jogosságának belátására csak az eloszlás fogalmához kell visszanyúlnunk:

E direkt formula használatára most csak két példát említünk.

Tegyük föl, hogy az X E[0,1) (azaz a [0,1) balról zárt jobbról nyílt intervallumban egyenletes) eloszlást követ, s szeretnénk az Y E[a,b)-belit segítségével előállítani. A direkt módszer szerint nincs más teendőnk, mint megkeresni a G(x)=P(Y<x) függvény inverzét. Tudva, hogy a G(x)=(x-a)/(b-a), az inverzére a G^-1(x)=x*(b-a)+a adódik.

Ugyanis G^-1(G(x))=G^-1( (x-a)/(b-a) )=x , amit teljesítenie kellett.

Második példaként az előbbi X-ből „származtassunk” egy Y m-paraméterű exponenciális eloszlású valószínűségi változót!

Az Y eloszlásfüggvénye: G(x)=1-e^-m*x.

Keressük azt a G^-1 függvényt, amelyre igaz a G^-1(G(x))=x ∀x-re, azaz G^-1(1-e^-m*x)=x.

Kérdés, hogy az 1-e^-m*x-ből milyen transzformációval „hozható ki” az x?

Ránézésre adódik:

Mivel 1-x eloszlása ugyanaz, mint az x-é, ezért igaz az (egyszerűbb) formula:

Ezen az ún. direkt módszeren kívül több más módszer is ismert. Egy-kettőre az alkalmazások során ki is térünk.

C. A kísérletezés módszertani alapja

1.-ben a folyamatot „állóképekre” bontottuk úgy, hogy a következő képet a folyamat változása determinálja. A 2.-ban megmutatott módon e változást valószerűvé tudtuk tenni, bármilyen eloszlást is kövessen az. Most megmutatjuk, hogy a szimulációs eljárásunk „filozófiája” – ti. hogy adott eloszlás szerint egyedi véletlen eseményeket generálunk, s ezek együttesen alakítják ki a valóságos képet – helyes. Ha ugyanis az adott eloszlás szerint generált eseményeket számláljuk, és képezzük ezek relatív gyakoriságát (pontosabban az empírikus eloszlást), akkor a kísérletek (esemény-generálások) számát növelve egyre jobb az egyezés az empírikus és az elméleti eloszlás között (ezt fejezi ki a Glivenko-tétel).

A tétel súlyát és „kényszerítő erejét” lássuk be oly módon, hogy számoljuk ki rögzített ε (ε>0) esetén mi a valószínűsége a

eseménynek egyre növekvő n-re. (Itt h(A,n) az A esemény gyakoriságát jelenti n kísérlet után, a P(A) az A valószínűsége.)

D. Meddig folytassuk a kísérletet?

A szimuláció során nem kísérletezhetünk végtelenségig. Egyrészt időkorlát miatt, másrészt a számítógép véges pontossága miatt sincs értelme bizonyos határon túl a kísérletek számát növelni. Van-e arra vonatkozó információnk, hogy meddig érdemes a kísérleteket végezni? A válasz igenlő, de ahogy már megszoktuk, valószínűségi jellegű. Ugyanis szimulációnál csak valamilyen valószínűséggel tudunk bármit is „garantálni”. Az előírt pontosságú közelítéshez szükséges kísérletek számát legegyszerűbben a Csebisev tétel alapján konstruálhatjuk meg:

Meghatározandó a kísérletek azon száma (n), amelyre a pontosság, 1-p a megbízhatóság:

A Csebisev tételt alkalmazva kapjuk:

Innen a feltétel kielégíthető, ha

Tegyük föl, hogy addig kísérletezünk, ameddig a kedvező esetek száma eléri az m-t: h(A,n)=m, ekkor a

érdekel bennünket!

Vissza a tartalomjegyzékhez

2.3. A szimuláció gyakorlata

véletlenszám generátorunk – feltevésünk szerint – [0,1) egyenletes eloszlást követ. Pontosabban a számítógép véges pontossága miatt ezt helyesebb lenne „diszkrét egyenletes” eloszlásnak nevezni, hiszen ha a pontossága d, akkor a véletlenszámok a következők lesznek: 0, d, 2*d, ..., i*d,..., N*d, ahol N:=1/d-1. Feltesszük, hogy a generátorunk a fenti számokat egyenletesen, azaz azonos d valószínűséggel állítja elő. Továbbá feltesszük, hogy egyéb számolásokra és konstansokra is jellemző a pontosság fenti d értéke.

Azt a kérdést kell megválaszolnunk, hogy milyen alapon várunk el a számítógépes szimulációtól a valóságra utaló, a szimulációs folyamatra jellemző adatokat, információkat, illetve hogyan kell egyáltalán komplex jelenségeket szimulálni. Keressük előbb az első kérdésre a választ!

Tekintsük a véletlenszám-generátort valószínűségi változónak (R), amely a valós [0,1) intervallumon van értelmezve és a d[0,1):={i*d: i=0..N} az értékkészlete.

A véletlenszám-generátor, mint függvény, diagramja.

Jelöljük R(i)-vel az i.-ként generált véletlenszámot. Ekkor az R(i) sorozat sztochasztikus folyamatot alkot.

Mekkora hibát követünk el akkor, ha eltekintünk a véletlenszám-generátor „diszkrét” voltától? A hibabecslés érdekében számítsuk ki e diszkrét egyenletes eloszlás várható értékét és szórását!

Közben felhasználtuk az N és a d közötti összefüggést, és a diszkrét egyenletességet: P(R = i * d) = d.

Azaz az egyenletestől való eltérés – a várható érték szempontjából – az ábrázolási pontosság felével egyenlő.

Hasonló következtetésre jutottunk a szórásnégyzetet illetően: a szórásnégyzet a pontosság négyzetének (!) 12-ed részével tér el a (folytonos) egyenletes eloszlású valószínűségi változó szórásnégyzetétől.

Tehát a levonható tanulság: Az eltérés a „valódi” egyenletes eloszlástól nem haladja meg az ábrázolási pontosság értékét, amiből látszik, hogy a bizonytalanság legfeljebb az utolsó jegyben mutatkozhat.

Összetett események szimulációja

A komplex jelenségek szimulációjának „hogyanját” közelítsük meg a legegyszerűbb alapeset felől! Egy p valószínűségű A eseménnyel kapcsolatos kísérletet kell lejátszanunk. Definiáljunk egy valószínűségi változót, X(A)-t, amely akkor 1, ha az A bekövetkezett, és 0 egyébként (X(A) az A esemény indikátorváltozója). Ennek tehát értelmezési tartománya az {A,-A} halmaz, és értékkészlete a {0,1}.

Jelölje R (valószínűségi változó) a véletlenszám-generátor adta számsorozat egy elemét! Miként kapcsolhatjuk össze X(A)-t és R-t úgy, hogy azt állíthassuk, hogy e két valószínűségi változó ekvivalens valószínűségi értelemben? A kapcsolatot egy függvény segítségével építjük ki. Első lépésként az R és az X(A) eltérő értelmezési tartományát kell megvizsgálnunk. Osszuk Do(R) pontjait két osztályba, s feleltessük meg az egyik osztályt A-nak, a másikat -A-nak. Definiálunk egy R(A) függvényt úgy, hogy Do(R(A))=Ra(R), Ra(R(A))= {0,1} és R(A)(R(o)):=1, ha az o az A-hoz lett rendelve, és 0, ha az o -A-hoz tartozik.

A valóságban elvégzett kísérletsorozatot jellemezzük az A bekövetkezés relatív gyakoriságának várható értékével. Hasonlóan a szimulációbeli kísérletben. Ha P(X(A)=1)=P(R(A)=1), akkor elfogadjuk ekvivalenseknek.

Itt kell folytatni tovább az R(A) valószínűségi változó értelmezését, azaz a Do(R) osztályfelbontását úgy, hogy P(R(A)=1) = P(A) teljesüljön.

Legyen

azaz Do(R) A-hoz rendelt osztálya!

Felhasználtuk, hogy P(R=i*d)=d.

Mivel a véletlenszám-generátorunk d-finomságú, így a véletlenszámot egyértelműen meghatározza az őt tartalmazó [i*d,(i+1)*d) részintervallum. Tehát vagy az adott osztályba kell tartozzon a teljes részintervallum, vagy diszjunkt tőle. Következésképpen P(R(A)=1|R=i*d)=1 vagy 0.

Itt az I azon indexek halmaza, amely részintervallumon a fenti valószínűség 1.

Vagyis olyan Do(R)-et kell kiválasztani, amely részintervallumainak száma:

Erre a legkézenfekvőbb a [0,P(A)) intervallumot kiválasztani (bár éppúgy megfelelő lenne az [1-P(A),1) választás is), hiszen részintervallumainak száma a fenti. (Megjegyezzük, hogy akár nem érintkező részintervallumok egyesítése által meghatározott halmaz is megfelelne, ha a fenti számú részintervallumból áll. Ez számítástechnikai szempontból azonban nem érdekes, hisz túl sok információt kell a halmazhoz tárolni.)

Összefoglalóan az R(A):[0,1){1,0},

Nos, az R(A) valószínűségi változó már magától értetődően előáll az R segítségével, csak az R(A)=1 és R<P(A) események azonosságát kell észrevennünk. (Bevezetjük a p:=P(A) jelölést.)

A. A és B diszjunkt események, együtt teljesen lefedik a biztos eseményt

Legyen az A esemény P, a B esemény pedig Q valószínűségű! Ha R<p, akkor az A eseményt bekövetkezettnek tekintjük, egyébként pedig a B-t.

Vegyük a [0,1) intervallumot, osszuk két részre:

Ha egy véletlenszám az első felébe esik, akkor az A esemény következik be, ha nem, akkor pedig a B.

Ritkán kell ilyen egyszerű esetet szimulálni, de sokszor ehhez hasonló módon, elvében megegyezően építhető föl a program szimulációs magja. Néhány ilyen esetet tekintünk át.

B. A és B nem diszjunkt események, de együtt teljesen „lefedik” a biztos eseményt

Az AB=0 nem teljesül, ezért nem alkalmazható az 1. módszer. Azonban könnyen megteremthető a kapcsolat 1-gyel, ha az {A,B} helyett az {A’:=A-AB, C:=AB és B’:=B-AB} eseményrendszert tekintjük. Az {A’,C,B’} már teljes eseményrendszert alkot.

Tehát

• A’ következik be, ha R<P(A’)=P(A-AB)=P(A)-P(AB),
• C,ha P(A’)≤R<P(A’)+P(C)=P(A),
• B’,ha P(A’)+P(C)≤R<P(A’)+P(C)+P(B’)=P(A)+P(B)-P(AB).

Mivel A+B=Ω (Ω a biztosan bekövetkező esemény), így P(A)+P(B)-P(AB)=1 egyenlőségnek igaznak kell lennie, következésképpen P(AB)=P(A)+P(B)-1.

Ezen eredményekből már könnyen származtathatók az alábbiak.

A,B: AB≠0, A+B=Ω események szimulálása:

• R<1-P(B) esetén az A egyedül következett be,
• 1-P(B) ≤R<P(A) esetén mindkét esemény,
• egyébként a B egyedül.

Jelölje P az A, Q pedig a B esemény valószínűségét! A fentiek szerint a [0,1) intervallumot az alábbiak szerint osztjuk három részre:

C. A és B diszjunktak, de nem adják ki a biztos eseményt

Az A+B=Ω feltétel nem teljesülése miatt megint nem alkalmazható az 1. Kiegészítve az {A,B} rendszert egy fiktív C (=Ω-A-B) eseménnyel: már készen vagyunk.

A,B: AB=0, A+B≠ Ω események szimulációja:

• R<P(A) esetén az A,
• P(A)≤R<P(A)+P(B) esetén a B következik be,
• egyébként egyik sem.

Jelölje P az A, Q pedig a B esemény valószínűségét! A fentiek szerint a [0,1) intervallumot az alábbiak szerint osztjuk három részre:

D. Az A, B és C eseményekről a következőket tudjuk:

A⊆B+C, és B, C események diszjunktak, és ismert a P(B), P(C) és P(AB). Ekkor elemi úton megkapható a P(A), a P(BA) és a P(B-A), ahol a -A az A komplemensét jelöli:

• P(A)=P(A(B+C))=P(AB)+P(AC) (mivel AB+C és BC=0)
• P(B|A) =P(AB)/P(A)
• P(B|-A)=P(-AB)/P(-A)=(P(B)-P(AB))/(1-P(A))(mivel AB⊆B).

Ekkor

• R<P(-AB) esetén csak a B következett be,
• P(-AB)≤R<P(B) esetén az A is és a B is,
• P(B) ≤R<P(B)+P(AC) esetén az A is és C is,
• P(B)+P(AC) ≤R<P(B)+P(C) esetén csak a C teljesül,
• egyébként egyik sem.

Az eddigi eseteket egyszerűségük és alapvető voltuk miatt nevezhetjük a szimuláció elemi eseteinek. A most következők már összetettebbek, de mindegyikük magától értetődően adódik az elemiekből.

E. n eseményből álló teljes eseményrendszer, azonos valószínűséggel

N eseményünk van, kizáróak, rajtuk kívül más nem lehet, egyenlő valószínűségűek. Ekkor minden esemény pontosan 1/N valószínűségű.

Ekkor az esemény indexe könnyebben meghatározható. Legyen P:=P(A_k)=1/n ∀k-ra!

Mivel

azaz

Az egyenlőtlenséget n-nel szorozva kapjuk, hogy

esetén kell bekövetkeznie az A_i eseménynek. Tehát az esemény sorszáma így határozható meg: i:=[n*R+1].

• 0≤véletlenszám<1
• 0≤N*véletlenszám<N
• 1≤N*véletlenszám+1<N+1
• 1≤egészrész(N*véletlenszám+1)≤N

F. n eseményből álló teljes eseményrendszer

A továbbiakban felhasználjuk a tárgyalt alapeset azon fontos eredményét, hogy a szimulálni kívánt eseményhez (ott az A-hoz, illetve a -A-hoz) a [0,1) intervallum megfelelő felosztását kell csak megkeresni.

Az A₁,...,A_n események teljes eseményrendszer volta azt jelenti, hogy

és i≠k esetén

Tehát minden esetben pontosan (!) 1 esemény következik be.

Ebből szemre is látszik, hogy olyan I_i részintervallum-felosztások között lelhető meg a keresett, amelyekre teljesülnek az alábbi feltételek:

és i≠k esetén

Tehát az első feltétel a b_i≥a_i+1-t kívánja meg, másrészt a második a b_i≤a_i+1-t vonja maga után, amik csak a b_i=a_i+1 esetben elégíthetők ki egyszerre.

Ez a részintervallum-partícionálás még csak a fenti két feltételt veszi figyelembe, semmi köze sincs az események valószínűségeihez. Az a_i-k megadása úgy kell történjen, hogy az [a_i,b_i) intervallum hossza P(A_i) legyen – ahogy ezt korábban beláttuk –, azaz b_i-a_i=P(A_i), s a teljesség miatt nyílván a₀=0, valamint b_n=1.

Következésképpen

Vegyük észre, hogy tulajdonképpen most is a korábban említett direkt módszert alkalmaztuk, csak épp egy „diszkrét függvény” inverzét kellett képeznünk! Nevezetesen az

függvényét.

Ahhoz, hogy a szimuláció eredményét helyesnek tarthassuk, az eseményeknek eleget kell tenniük az alábbi két követelménynek:

• ∀A_i: P(A_i)>d (a P(A_i)=0 eseményt a teljes eseményrendszerből kizárjuk, mint érdektelent),
• n≤1/d. Az egyenlőség is csak az 1/n elemi valószínűségű egyenletes eloszlás esetén lehetséges. (Ez utóbbi feltétel lényegesen befolyásolja a „folytonos” szimulációt!)

Az n-elemű teljes eseményrendszer szimulációja történhet így: a szimuláció során az A esemény következett be, ha

Utaljunk ennek „praktikus” (=program) megvalósítására is! A megfelelő esemény kiválasztásakor nem célszerű az eseményazonosítás minden egyes alkalmával a feltételben szereplő két szummát újra meg újra kiszámolni, hanem járjunk el így:

Az algoritmusbeli P vektor az A_i események valószínűségeit tartalmazza.

G. Végtelen sok esemény

Végtelen sok eseményünk van, kizáróak, rajtuk kívül más nem lehet, különböző (P_i) valószínűségűek. Ekkor adott tetszőleges ε>0 számhoz legyen m a legnagyobb index, amire P_m>ε!

H. Binomiális eloszlás

Az n független kísérlet során az A esemény bekövetkezésének gyakoriságait {h(A,n)=k} tekintsük elemi eseményeknek! Állítsuk elő ezeket az elemi eseményeket! (A p:=P(A) és n paraméterű binomiális eloszlás generálása.)

Ekkor az A esemény bekövetkezéseinek gyakorisága határozza meg az elemi események halmazát: A_i:={h(A,n)=i}.

A hozzájuk tartozó valószínűségek:

Az A_i teljes eseményrendszer, így az F.-re való hivatkozással elintézettnek is tekinthetjük. Mégis, mielőtt lezárnánk e kérdést vegyük szemügyre az F. alatt szereplő megjegyzéseket!

• P(A_i)≥d, itt ezzel azonos feltétel a min(P(A₀),P(A_n))≥d (hiszen P(A_i)>min(P(A₀),P(A_n)) ∀i-re). Az általánosságot nem sértjük meg azzal, hogy feltesszük: p≤1/2. Ekkor min(P(A₀),P(A_n))=P(A_n)=p_n, tehát p_n≥d kell igaz legyen, tehát a még használható n-re kapjuk az ln(d)/ln(p) felső korlátot.
• az n≤ln(d)/ln(p)-t összehasonlítva az n-re kapott korábbi n≤1/d-vel láthatóan lecsökkentette a lehetséges n-ek körét. (Ugyanis 1/d>>ln(1/d)/ln(1/p), ha p<<1).

Az első eset programozásához a következőket érdemes észrevenni:

• Ha a h(A,n) binomiális eloszlású valószínűségi változó szimulálását az F. alatti módon végezzük, akkor a valószínűségi változó egyes megvalósulásaihoz tartozó valószínűségeket csak egyszer célszerű kiszámolni és tárolni egy vektorban – mégpedig a valószínűségi változót megvalósító szubrutin első hívását megelőzően –, s később már csak a teljes esemény rendszerből kell kiválasztani egy elemét, melynek indexe lesz épp a h(A,n) „értéke”.
• Egy másik módszer lényege az, hogy az esemény indexe n db E[0,p)-beli véletlenszám összegének egészrészeként állítható elő.

Binomiális eloszlású véletlenszám előállítható a definíció szó szerinti alkalmazásával is – hányszor következik be egy p valószínűségű esemény n kísérletből:

I. Geometriai eloszlás

Egy p valószínűségű esemény első bekövetkezésének sorszámára vagyunk kíváncsiak:

Láthatólag alkalmazható a G. eset, van végtelen sok, egymást kizáró, különböző valószínűségű eseményünk.

Itt is egyszerűbb azonban a definíció szó szerinti alkalmazása:

Belátható,hogy egy véletlenszámra a

összefüggés éppen p(n) valószínűséggel igaz, azaz

J. n állapotú rendszer szimulációja

Legyen az R rendszer állapotainak halmaza az S véges halmaz. A rendszer adott s_i∈S állapotból valamely s_k∈S állapotba p_ik valószínűséggel halad. Írjuk meg a rendszer állapotát szimuláló program magját!

Kézenfekvőnek látszik ezt is visszavezetni az A. esetre. Ez könnyen megvalósítható, ha az ott szereplő vektor helyett most mátrixot alkalmazunk az alábbi értelmezéssel: P:=[p_ik], p_ik:=P(k. állapot|i. állapot). A mátrix i. sora így egy adott állapotból való átmenetekhez rendeli a valószínűségeket.

Az n állapotú rendszer szimulációja történhet így, ha

akkor a szimuláció során a rendszer az i. állapotba került (feltéve, hogy ezt megelőzően az m.-ben volt).

Egy kis „elméleti” jellegű kitérő, gyakorlati következményekkel:

Gyakran előfordul, hogy a modellezettről még sztochasztikusan sem áll rendelkezésünkre elegendő információ. Néha ilyenkor is kínál megoldást a szimulációs megközelítés. Ha bizonyos feltételek teljesülése garantálható, akkor két ismert típusú eloszlásra vezethető vissza a probléma.

K. Poisson-folyamatként írható le a jelenség, ha

• Homogén növekményű (stacionárius), azaz adott időtartam alatt a jelenség bekövetkezés gyakoriságának eloszlása csak az időtartam nagyságától függ, azaz h(t,t+d) és h(s,s+d) eloszlása megegyezik. (A h(t₁,t₂):= a jelenség bekövetkezésének gyakorisága t₁ és t₂ időpontok között.)
• Független növekményű, utóhatásmentes, azaz a jelenség bekövetkezéseinek gyakoriságát nem befolyásolja, hogy addig hányszor következett be: h(t₁,s₁), h(t₂,s₂), ..., h(t_n,s_n) valószínűségi változók teljesen függetlenek, ha a [t_i,s_i) intervallumok diszjunktak.
• „Ritka”, azaz annak a valószínűsége elhanyagolható, hogy egy kellően rövid időtartam alatt egynél többször következzék be az esemény:
  
  
• Néha kiegészítik a fenti feltételeket egy továbbival is, ami a precíz levezetéshez valóban fontos, ám gyakorlati vizsgálódásaink során ezt természetessége miatt figyelmen kívül hagyhatjuk (ti. a jelenség nem elfajuló). Ez azt jelenti, hogy 0<P(h(t,t+d)=0)<1. (Érdemes meggondolni, hogy e feltétel, milyen szélsőséges eseteket zár ki!)

A feltételeket teljesítő jelenségek Poisson-eloszlással írhatók le (k=0,1,...,):

ahol λ a jelenségre jellemző pozitív valós szám.

A Poisson-eloszlás nagy jelentőségű voltát bizonyítja az is, hogy bizonyos feltételek mellett a binomiális eloszlás határértékeként is megkapható! Ez a tény annál is fontosabb, mivel – különösen nagy n-ek-re – a binomiális eloszlású valószínűségi változó szimulálása gyakorlati akadályokba (végrehajtási idő-, vagy helyproblémába) ütközik. Ilyenkor komoly gyakorlati haszonnal jár a Poisson-eloszlásra való át-térés lehetősége.

Állítás: rögzített k-ra, ha n→+∞ és n*p→λ>0, akkor

Hogyan lehet Poisson-eloszlású valószínűségi változót szimulálni?

Ahogy korábban láttuk, a megoldáshoz úgy is eljuthatunk – kiindulva az E[0,1) eloszlású véletlenszám-generátorból –, hogy képezzük a Poisson-eloszlás „inverz-függvényét”. A teljes eseményrendszernél alkalmazott ötletből kiindulva azt az n egész értéket kell megtalálnunk, amelyre igaz

A Poisson-eloszlású véletlenszámokat generáló függvény megírásánál az alábbi rekurzív képletből lehet kiindulni:

T(n):=T(n-1)*λ/n; S(n):=S(n-1)+T(n)

L. Gauss-folyamathoz való konvergencia

Független valószínűségi változók összege normális (Gauss-) eloszlást követ, ha a tagok száma minden határon túl növekszik. A pontosabb határeloszlás-tétel a következő:

Ha X₁,… X_n valószínűségi változók függetlenek és egyforma eloszlásúak, létező m várható értékkel és s szórással, akkor az

összeg aszimptotikusan standard normális eloszlású, vagyis

A véletlenszám-generátorunk esetén – mint láttuk – m=1/2 és s=1/gyök(12). Most már minden ismeret birtokunkban van ahhoz, hogy megadjuk azt az algoritmust, amellyel közelítőleg N(μ,σ) eloszlású véletlenszámot generálhatunk.

Egy G N(μ,σ) normális eloszlású véletlenszám előállítása:

Most rendelkezünk egy közelítőleg

eloszlású véletlenszámmal, már csak a kérdéses normális eloszlásúvá transzformálása maradt hátra:

Ha n=12*m², akkor:

Igazolható, hogy viszonylag kicsi n esetére (n<10) is elég jó a közelítés. A pontosság az összeadandók számának növelésével tovább javítható, persze a műveletigény rovására. Speciális módszerekkel az n növelése nélkül is elérhető javulás.

Ahogy a Poisson-eloszlásnál örülni lehetett annak, hogy helyettesíteni – pontosabban: közelíteni – tudtuk vele a binomiális eloszlást, úgy ez az öröm még teljesebb lehet a normális eloszlású véletlenszám szimulálására adott fenti egyszerű algoritmus és az alábbi tétel ismeretében.

Moivre-Laplace tétel:

Rögzített z₀ és z₁ valós számok esetén, ha n→+∞, akkor

E tétel gyakorlati haszna akkor nyilvánvaló, ha a binomiális eloszlás realizálása számolásigényessége miatt nem kivitelezhető, az egyszerűen megvalósítható normális eloszlásúval való közelítéssel szemben.

Vissza a tartalomjegyzékhez