Iránymenti integrálás – ortonormált becslések

integrál kiszámítását, ahol f(z) azX halmaz indikátorfüggvénye, azaz f(z) =

( 1, haz∈X, 0, egyébként.

A (3.4) egyenlőség jobboldala alapján a következő Monte Carlo eljárást lehet megadni a P r{X}valószínűség kiszámítására. Legyen aξ valószínűségi változóϕ sűrűségfüggvényű, és jelölje ennek a független mintáitx_i, i= 1, . . . , N. Ekkor I egy torzítatlan becslése a

átlag. Ezt a becslést a durva becslésnek, vagy elfogadás-elvetés becslésnek nevezzük, mivel a gyakorlatban ez egy olyan eljárásra vezet, ahol a x_i, i = 1, . . . , N vektorok generálása után már csak azX halmazban fekvő vektorok számát kell meghatározni, így a becslés a relatív gyakoriság.

Közismert, hogy egy normális eloszlású ξ vektor felírható

ξ =χ_nTη (3.6)

alakban, ahol χ_n egy n szabadságfokú χ-eloszlású skalár valószínűségi változó, T egy felső–háromszög alakú mátrix, amelyre T T⁰ = R fennáll, továbbá az η valószínűségi vektorváltozó egyenletes eloszlású az egységgömb felületén, az S = {x | Pn

i=1 x²_i = 1}

halmazon. A ξ vektor χ_n „hossza” és az η „iránya” független valószínűségi változók.

Jelölje a χ-eloszlású χ_n valószínűségi változó eloszlásfüggvényét K(λ), λ ≥ 0 és az η valószínűségi vektorváltozó eloszlásfüggvényét V(y),y ∈ S. Ezekkel a jelölésekkel az integrálunk átírható: Vezessük be a g(y) jelölést a belső integrálra, legyen

g(y) = Z ∞

0

f(λTy)dK(λ). (3.8)

Ez a g(y) függvény egy rögzített y esetén megadja a λTy, λ ≥0 sugár X halmazba eső részének a valószínűségi tartalmát. Tegyük fel most, hogy a λz =λTy egyenes elmetszi

azX konvex halmazt, és definiáljuk az egyenes belépési, illetőleg kilépési pontjait aλ_Lés K(λ_U | y) jelölésekkel is fogunk élni. (Jegyezzük meg, hogy az egydimenziós K(·) elos-zlásfüggvény könnyen kiszámítható létező szoftverek segítségével.) Bevezetve a λ⁺_L = max{0, λ_L}, λ⁺_U = max{0, λ_U} jelöléseket, látható, hogy a g(·) függvény a

g(y) =K(λ⁺_U|y)−K(λ⁺_L|y). (3.10) formába írható. Hasonlóan kapható a λ⁻_U = min{0, λ_U}, λ⁻_L = min{0, λ_L} jelölések használatával

g(−y) = −K(−λ⁻_U| −y) +K(−λ⁻_L| −y),

vagyis a λL, λU állandók meghatározása után aλz=λTyegyenes valószínűségi tartalma is meghatározható:

e(z) =e(Ty) = [g(y) +g(−y)]/2.

A belső integrálra kapott (3.10) kifejezést visszahelyettesítve a (3.7) kettős integrálba írhatjuk, hogy

Az egyenlet alapján a következő Monte Carlo eljárás adható azIkiszámítására: generáljuk az y_i, i = 1, . . . , N független mintákat az S-en adott egyenletes eloszlásból és számítsuk ki determinisztikusan (a megfelelő program segítségével) a g értékét minden y_i esetén.

Tehát becslésünk I-re a következő: Ez a becslés az I-t a következő módon határozza meg: a (3.8) belső integrált deter-minisztikusan számítjuk, míg a külső integrálra mintavételt hajtunk végre. A becslést hatékonyabbá tudjuk tenni a következő két módosítás segítségével.

Az egyik módosítást arra alapozzuk, hogy a Θ₂ szórása viszonylag nagy, mert az y_i vektorok „túl véletlenszerűen” vannak szétszórva az egységgömb felületén, így egy olyan elrendezésre van szükségünk, amely a felhasznált vektorok „egyenletességét” növeli. Ez a következő módon érhető el: a független yi vektorok helyett egy ortonormalizált vek-torokból álló rendszert állítunk elő.

Tekintsük a véletlen U rendszert, amely S-ben lévő, ortonormalizált vektorokból áll;

legyen U ={uⁱ, i= 1, . . . , n|uⁱ ∈S,uⁱ⁰u^j =δ_ij, i, j = 1, . . . , n}, ahol δ_ij = 0, i6=j, δ_ii=

1, és U egyenletes eloszlású az ortonormalizált rendszerek felett. A második módosítás előkészítéseként tekintsük tetszőleges két U-beli vektor összegét és különbségét, vagyis legyen

v^i,j,s = 1

√2 s₁uⁱ+s₂u^j

, (3.13)

ahol azi, j indexpár és az s előjelvektor az összes lehetséges értéket felveszi a

J^∗ ={(i, j,s)|i= 1, . . . , n−1, j = 2, . . . , n, i < j, s₁ = 1, s₂ = 1,vagy s₁ = 1, s₂ =−1}

halmazból. Az uⁱ vektorok normalizált összegét azért használjuk, hogy az egy vek-torra eső számítási munkát csökkentsük. Egy adott U rendszer esetén csak n számú uⁱ vektorunk van, de az U-ból előállítható v^i,j,s vektorok száma 2n(n −1) (az előállított egyenesek száma pedign(n−1)).

A másik módosítást a következő, számítástechnikai szempontból fontos megjegyzésre alapozzuk. A Tv^i,j,s transzformált vektorok kiszámítása helyett csak azuⁱ vektorok Tuⁱ transzformált vektorait határozzuk meg, mivel ezek segítségével az előbbiek megadhatók, ugyanis

z^i,j,s =Tv^i,j,s = 1

√2 s₁Tuⁱ +s₂Tu^j ,

tehát n(n − 1) mátrixszorzás helyett csak n mátrixszorzásra lesz szükségünk a Tv^i,j,s vektorok előállításához. Így az egy generált vektorra eső számítási munka lényegesen csökkenthető, továbbá a generált vektorok „egyenletessége” is növelhető.

A két módosítás figyelembevételével megadjuk egy teljes U rendszer esetén a becslés formális alakját akkor, ha tetszőleges kétU-beli vektor összegére és különbségére végezzük el a függvényértékek kiszámítását: Ezt a becslést ortonormalizált–2 (ortonormált–2), vagy röviden O₂ becslésnek nevez-zük. Természetesen a gyakorlati számítások során N különböző U rendszert generálunk véletlenszerűen, és az eredményül kapott valószínűségek átlagát számítjuk ki.

Két darab, U rendszerből származó vektor összege és különbsége helyett vehetnénk k darab U-beli vektor minden lehetséges előjellel vett normalizált összegét is, ezáltal az O_k becslést kapnánk – ezzel megnövelnénk az egy rendszerből előállítható vektorok számát és csökkentenénk az egy előállított vektorra eső számítási munkát. Az így kapott Ok becslések formális leírása könnyen megkapható a fentiek alapján, így ezek leírását mellőzzük. Az O₂ becslések általában számítástechnikailag megfelelőnek bizonyultak (az O_k becslések összehasonlítására vonatkozó numerikus eredmények találhatók például a [De 79], [De 98c] cikkekben).

Az O₂ becslés végleges leírásában a λ⁺_L, λ⁺_U és a g függvény helyett az eredeti λ_L, λ_U állandókat és az e függvényt használjuk, vagyis egyidejűleg számítjuk a g(y) és g(−y) függvényértékeket – természetesen ez megfelezi azon vektorok számát, amelyekre aλ_L, λ_U értékeket ki kell számítani. Az O2 becslést megvalósító algoritmus lényegi lépéseit az alábbiakban adjuk meg.

O₂ algoritmus (általános eset, egy U rendszer esetén)

1. Generáljuk az U ={uⁱ} rendszert és számítsuk ki az uⁱ =Tuⁱ, i= 1, . . . , n vektorokat.

2. Legyen Sum= 0.

3. Határozzuk meg az összes lehetséges z=z^i,j,s=Tv^i,j,s = ^√1

2 s₁uⁱ+s₂u^j , (i, j,s)∈J^∗ vektorokat és minden z vektorra végezzük el a 4. lépést.

4. Kezdet: (az e(z) függvény kiszámítása)

kiszámítjuk a λ_L = min{λ |λz∈X}, λ_U = max{λ|λz∈X} értékeket, ha nincs metszés, akkor legyen λ_L=λ_U = 0,

ha λU ≥λL ≥0, akkor legyen Sum=Sum+K(λU)−K(λL), ha λ_U ≥0≥λ_L, akkor legyen Sum=Sum+K(λ_U) +K(−λ_L), ha 0≥λ_U ≥λ_L, akkor legyen Sum=Sum+K(−λ_L)−K(−λ_U), vége (az e(z) függvény kiszámítása).

5. Adjuk át a keresett valószínűség torzítatlan becsléseként a p=Sum/[2n(n−1)] értéket.

Megjegyezzük, hogy egy teljes algoritmusban az U rendszernek N darab független re-alizációját használjuk, N-et nevezzük a mintaszámnak. Az O₂ becslés kiszámítása azért hatékonyabb a durva módszernél, mert egyU rendszer eseténO(n²)mátrixszorzás helyett csakO(n)mátrixszorzásra van szükség és még néhány skalár szorzásra. Ez a számítástech-nikai hatékonyság növekedés az összes, továbbiakban leírásra kerülő esetben igaz, vagyis eloszlásfüggvény, téglatest, poliéder, ellipszoid, körkúp valószínűségének kiszámítása es-etén is.

A (3.7) egyenletben leírt integráltranszformációt és az eredményül kapott kettős inte-grált többféleképpen is felfoghatjuk.

(i) Tetszőleges p=P{ξ ∈X}valószínűség felírható, mint az X halmaz f(ξ) =

( 1, haξ ∈X,

0, egyébként (3.15)

indikátor valószínűségi változójának a várható értéke, vagyis p=E[f(ξ)].

Ez a várható érték megfelel a (3.4) jobboldalán álló integrálnak. Felhasználva aE(α) = E[E(α|β)] ismételt (feltételes) várható érték összefüggést ez a várható érték átírható a

p=E[f(ξ)] =E[f(χ_nTη)] = E[E(f(χ_nTη)|η)] (3.16)

alakba, ami viszont pontosan megfelel (3.7) kettős integráljának.

(ii) Egy másik lehetséges értelmezés adódik a numerikus integrálás szempontjainak figyelembevételével. A Monte Carlo integrálás elég jól működik, ha a feladat dimenziója nagy, de tudjuk, hogy viszonylag lassú, O(N^−1/2) a konvergencia sebessége. A hagy-ományos (determinisztikus) integrálási szabályok kis hibával dolgoznak, de ezeket nem nagyon lehet magasabb dimenzióban használni.

A kettős integrál formájába írt kifejezés a munkánkat két részre osztja: egy egydi-menziós, vonal menti integrál meghatározása hagyományos numerikus integrálási tech-nika segítségével és egy, az n-dimenziós térben elhelyezkedő (n−1)-dimenziós felületen elvégzett Monte Carlo integrálásra. Ezt a felbontást sugaras-felületi (radial–spherical) integrálásnak [MG 97], vagy iránymenti szimulációnak (directional simulation) is nevezik [DB 89], s a többdimenziós t-eloszlás eloszlásfüggvényének kiszámítására, illetőleg más halmazok valószínűségének meghatározására is használható (elliptikusan szimmetrikus sűrűségfüggvények esetén).

(iii) Végül a Monte Carlo integrálás szempontjából is megvizsgálhatjuk a dekompozí-ciót. Minden szimuláció esetén a fő kérdés az, hogyan lehet csökkenteni a becslés szórását (anélkül, hogy lényegesen megnövelnénk a szükséges munkát). Ez az eljárás éppen erre példa – a változók számának csökkentésével szóráscsökkenést érünk el. A Monte Carlo in-tegrálás területén szokásos szóhasználattal egy ortonormált becslést egy determinisztikus integrálási formula randomizált változatának is nevezhetünk.