Eötvös Loránd Tudományegyetem Természettudományi Kar
MCMC módszerek Szakdolgozat
Sztahó Roland István Matematika B.Sc., elemz® szakirány
Témavezet®:
Pröhle Tamás,
tanársegéd
Valószín¶ségelméleti és Statisztika Tanszék
Budapest 2012
Tartalomjegyzék
1. Bevezetés
3
1.1.
Integrálás Monte Carlo módszerrel
. . . . . . . . . . . . . . . . . . . .
3
1.2.
Markov láncok . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
4
2. MCMC szimuláció
8
3. Metropolis-Hastings módszer
10
3.1.
A Metropolis-Hastings algoritmus . . . . . . . . . . . . . . . . . . . . .
10
3.2.
Megvalósítás . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
12
4. Gibbs Sampling
14
4.1.
A Gibbs mintavételez® eljárás
. . . . . . . . . . . . . . . . . . . . . . .
14
4.2.
Konvergencia
. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
15
4.3.
Két- és többváltozós eset . . . . . . . . . . . . . . . . . . . . . . . . . .
17
4.4.
Megvalósítás . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
19
5. Slice sampling
20
5.1.
Egyváltozós eset . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
21
5.2.
Többváltozós eset . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
26
5.3.
Egyenletes eloszlás egy kétdimenziós tartományon . . . . . . . . . . . .
27
6. Perfect sampling
29
6.1.
Az alapötlet . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
29
6.2.
A Coupling From The Past (CFTP) algoritmus
. . . . . . . . . . . . .
30
6.3.
Fontainebleau lehulló levelei
. . . . . . . . . . . . . . . . . . . . . . . .
31
6.4.
Sandwiching . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
32
1
7. Összefoglalás
38
8. Függelék
39
8.1.
R program a sandwiching módszer bemutatásához . . . . . . . . . . . .
9. Köszönetnyilvánítás
39
45
2
1. fejezet Bevezetés
MCMC (Markov chain Monte Carlo) módszerek alatt algoritmusok egy családját értjük, melyek a második világháború alatt, Los Alamosban láttak napvilágot nem sokkal az eredeti Monte Carlo módszerek születése után. Nicholas Metropolis több társszerz®vel 1953-ban publikálta azt a módszert, melyet utána Metropolis algoritmusnak nevezünk. Ezt kés®bb W. Keith Hastings 1970-ben általánosította amelyet azóta is mint az egyik leismertebb eljárást, Metropolis-Hastings algoritmusnak nevezzük. A Geman testvérek 1984-ben publikált cikkükben bemutatták a Metropolis-Hastings algoritmus egy speciális esetét, amely Gibbs mintavételez® eljárás néven híresült el. A számítógépek fejl®désével a módszereket egyre szélesebb körben alkalmazták, de a statisztikában csak 1990 után váltak igazán népszer¶vé.
Az algoritmus-családnak számos alkalma-
zási területe van, például: numerikus integrálás, Bayes-i statisztika, statisztikus zika, bioinformatika.
1.1.
Integrálás Monte Carlo módszerrel
Az eredeti Monte Carlo módszert zikusok fejlesztették ki komplex integrálok kiszámítására véletlen számok segítségével. Tegyük fel, hogy szeretnénk kiszámítani a következ® integrált:
Z
b
f (x)dx. a Tegyük fel, hogy ahol
h(x)
egy
(a, b)
f (x)-et
fel tudjuk írni a következ® alakban:
intervallumon értelmezett s¶r¶ségfüggvény.
3
f (x) = g(x) · h(x),
Ekkor
Z
b
b
Z
g(x) · h(x)dx = Eh(x) [g(x)].
f (x)dx = a
a
Tehát az integrált kifejezhetjük mint
g(x). n
Ezáltal ha van egy
x1 , . . . , x n
h(x)
mintánk a
várható értékét, ha a s¶r¶ségfüggvény
g(x)
s¶r¶ségfüggvényb®l, kell®en nagy
esetén
Z a
b
n
1X f (x)dx = Eh(x) [g(x)] ' f (xi ). n i=1
Ezt a közelítést nevezzük Monte Carlo integrálásnak.
1.2.
Markov láncok
Miel®tt rátérnénk különböz® MCMC módszerek tárgyalására, a Markov láncokkal kapcsolatos néhány alapfogalmat elengedhetetlen tisztázni. A Markov láncok szemléletes értelmezéséhez vegyünk egy véges sok csúcsból álló, élsúlyozott, irányított gráfot, melyben hurokélek is megengedettek. Legyen minden él súlya nemnegatív, valamint minden csúcs esetén a kimen® élek súlyainak összege legyen
1.
Ha ezen a gráfon bolyongunk
az élekre írt valószín¶ségek szerint, valamint egységnyi id®közönként lépünk, akkor diszkrét paraméter¶ Markov láncot kapunk.
1.1. De níció. Legyen adott az {Xn }n∈N folyamat, ahol Xn : Ω → I valószín¶ségi
változók az (Ω, F, P ) valószín¶ségi mez®n, I pedig megszámlálható halmaz. A folyamatot diszkrét paraméter¶ homogén Markov láncnak nevezzük, ha a folyamat Markov tulajdonságú, azaz Fn = σ(X0 , X1 , . . . , Xn )
jelöléssel minden B ⊆ I és minden m ≥ n esetén P (Xm ∈ B | Fn ) = P (Xm ∈ B | Xn ).
A Markov folyamat stacionárius (vagy homogén) átmenetvalószín¶ség¶, azaz minden i, j ∈ I esetén minden olyan n-re, melyre P (Xn = i) > 0, P (Xn+1 = j | Xn = i) = pij ,
n-t®l függetlenül. 4
Az I halmazt, azaz a valószín¶ségi változók lehetséges értékeit állapottérnek nevezzük. Tehát a Markov lánc egy olyan Markov tulajdonságú, azaz a
X0 , X1 , . . . , Xn
(n + 1)-edik
valószín¶ségi változó sorozat, mely
pillanatbeli állapot meghatározásához az
el®z® állapotok nem szükségesek, csupán az éppen aktuális állapotok ismerete nem változtatja meg egyik A
P = (pij )
X0
voltunk.
ugyanis a korábbi
átmenet valószín¶ségét sem.
mátrixot átmenetmátrixnak nevezzük, elemei az átmenetvalószín¶sé-
gek, azaz hogy az
i-ben
i→j
n-edik,
(n + 1)-edik
id®pontban
j -ben
vagyunk, feltéve hogy az
n-edikben
Az átmenetvalószín¶ségek összességét átmenetmagnak nevezzük.
eloszlását kezdeti eloszlásnak hívjuk, és
p = (pi )-vel
Az
jelöljük. Ez a két objektum
már meghatározza a folyamatot, hiszen a véges dimenziós eloszlásokat és a Markovtulajdonságot felhasználva:
P (X0 = i0 , X1 = i1 , . . . , Xn = in ) = = P (X0 = i0 )P (x1 = i1 | X0 = i0 ) . . . P (Xn = in | X0 = i0 , . . . , Xn−1 = in−1 ) = = pi0 pi0 i1 . . . pin−1 in . Jelölje
πj (t) = P (Xt = j) annak a valószín¶ségét, hogy a lánc a a sorvektort, melynek
i-edik
eleme
melynek gyakran minden eleme ladtával ez az
1
0,
t-edik
πi (t).
id®pontban
A láncot egy
j -ben
π(0)
van és jelölje
π(t)
azt
kezd®vektorral indítjuk,
egy kivételével, mely ebb®l adódóan
valószín¶ség eloszlik
π(t)
1.
A
t
el®reha-
koordinátáin.
1.2. De níció. A P mátrix
1. sztochasztikus, ha pij ≥ 0 és minden sorban az elemek összege 1, 2. duplán sztochasztikus, ha sztochasztikus és minden oszlopban az elemek összege 1, 3. szubsztochasztikus, ha pij ≥ 0 és minden sorban az elemek összege legfeljebb 1. 1.1. Tétel. Tetsz®leges I megszámlálható halmazon adott p eloszláshoz és | I | × | I | méret¶ P sztochasztikus mátrixhoz létezik I állapotter¶ Markov lánc, melynek kezdeti eloszlása p, átmenetmátrixa P . 5
Annak a valószín¶sége, hogy egy lánc a
(t + 1)-edik id®pillanatban az si
állapotban
van, meghatározható a Chapman-Kolomogrov egyenl®ség alapján:
πj (t + 1) = P (Xt+1 = si ) = X
P (Xt+1 = si | Xt = sk ) · P (Xt = sk ) =
k
X
pki · πk (t).
k
Ezt könnyen láthatjuk, ha iteratívan alkalmazzuk az egyenl®séget:
π(t) = π(t − 1)P = (π(t − 2)P )P = . . . = π(0)P t . Ezek alapján de niáljuk az
n-edrend¶
átmenetvalószín¶séget:
1.1. Állítás. Legyenek pnij = P (Xn+m = j | Xm = i) az n-edrend¶ átmenetvalószí-
n¶ségek (feltesszük, hogy P (Xm egyenl®ség:
= i) > 0).
pn+m = ij
Ez azt jelenti, hogy a
X
Ezekre teljesül a Chapman-Kolomogrov
(n) (m)
pik pkj .
k
pnij
mennyiségek éppen a P n mátrix megfelel® elemei.
Szükségünk van még továbbá az állapotok osztályozására:
1. Azt mondjuk, hogy az i állapotból elérhet® j , (i → j ), ha van olyan n ≥ 0, hogy pnij ≥ 0.
1.3. De níció.
2. Azt mondjuk, hogy i és j közlekednek (lehet hogy több lépés szükséges), ha i-b®l elérhet® j és fordítva. Ez ekvivalenciareláció, tehát osztályokra bontja az állapotteret (csak az átmenetmátrixtól függ). A Markov lánc irreducibilis ha csak egyetlen osztályból áll. Ezek a tulajdonságok osztálytulajdonságok, azaz vagy egy osztály minden eleme rendelkezik az adott tulajdonsággal, vagy egyik sem.
1.4. De níció. Az {n > 0 : p(n) ii > 0} halmaz legnagyobb közös osztója az i periódusa,
jele: d(i). Ha a halmaz üres, akkor a periódust nem értelmezzük. Ha d(i) = 1, akkor az állapotot aperiodikusnak nevezzük. 6
Szemléletesen ez azt jelenti, hogy semmilyen x lépésszám esetén nem feltételezhetünk törvényszer¶ visszatérést egy adott állapotba.
1.5. De níció. Legyen P egy átmenetmátrix. A (pi )i∈I eloszlás stacionárius, ha pi =
minden i ∈ I esetén, azaz a pi kezdeti eloszlású, P átmenetvalószín¶ség¶ Xn Markov láncra P (Xn = i) = pi minden i ∈ I, n ≥ 0 esetén. Ez utóbbi esetben azt mondjuk, hogy a Markov lánc stacionárius. P
pk pki
k∈I
Ebb®l adódóan, ha
π ∗ -ra,
mint a
P
π∗
stacionárius eloszlás, akkor
mátrix baloldali
λ=1
π∗ = π∗P .
Így tekinthetünk
sajátértékhez tartozó sajátvektorára. A sta-
cionárius elolszás feltétele, hogy a Markov láncunk irreducibilis és aperiodikus legyen. Ha egy Markov lánc periodikus, akkor megvan annak a lehet®sége, hogy két állapot között váltakozik és nem tart egy stacionárius eloszláshoz, illetve a váltakozás maga lesz a stacionárius eloszlás. Ha
P -nek nincs (−1)-gyel egyenl® sajátértéke, akkor a lánc
aperiodikus. Az egyértelm¶séghez szükség van továbbá a megfordíthatósági feltételre:
pij πi∗ = pji πj∗
∀i, j.
1.6. De níció. A fenti egyenl®ség teljesülése esetén a Markov láncot megfordíthatónak
nevezzük.
Könnyen láthatjuk, hogy ebb®l már következik
(πP )j =
X
πi pij =
X
i
πj pji = πj
π = πP ,
X
i
ugyanis
πP j -edik
eleme:
pji = πj · 1 = πj .
i
A diszkrét állapotter¶ Markov láncok általánosítása folytonos állapottérre azon a megfontoláson alapszik, hogy az átmenetmag,
P (x, y)
kielégíti az
Z P (x, y)dy = 1 egyenl®séget, illetve a Chapman-Komologrov egyenl®séget értelmezzük folytonos esetben:
Z πt (y) =
πt−1 (x)P (x, y)dy.
Így (ha létezik) a stacionárius eloszlás, kielégíti az alábbi egyenl®séget:
∗
π (y) =
Z
π ∗ (x)P (x, y)dy.
7
2. fejezet MCMC szimuláció
Ahogy láttuk az el®z® fejezetben, a Markov láncok általánosságban a követkz® módon m¶ködnek: adott egy átmenetmag, mely tartalmazza a (feltételes) átmenetvalószín¶ségeket arra vonatkozóan, hogy egy pontból egy másikba lépünk, legyen ez ahol az
x ∈ Rd , A
pedig
σ -algebra.
P (x, A),
Az átmenetmag egy feltételes eloszlásfüggvény, mely
x állapotból egy másik állapotba lépés valószín¶ségét írja le A-n, ezért P (x, Rd ) = 1,
illetve mivel megengedjük hogy az
x állapotból x-be lépjünk, P (x, {x}) nem feltétlenül
nulla. A célunk az volt, hogy meghatározzuk azokat a feltételeket melyek esetén létezik a
π ∗ (x)
stacionárius eloszlás, melyhez az iteráció során tart az átmenetmag, amely
kielégíti a következ® egyenl®séget:
Z π(y) =
P (x, y)π(x)dx. Rd
Az MCMC módszerek pont fordítva m¶ködnek: a kívánt az átmenetmag ismeretlen. Ahhoz, hogy
∗
π (·)-ból
π ∗ (·)
eloszlás ismert, de
mintákat generáljunk, olyan átme-
netmagot kell gyártanunk, mely az iterációk számának növelésével konvergál a
π ∗ (·)
ismert eloszláshoz. A folyamatot egy tetsz®legesen választott vektorral kezdjük és kell®en sok iteráció után a szimuláció során meg gyelt értékek eloszlása közelít®leg azonos lesz a
π ∗ (·) ismert eloszlással.
Ez nagy segítséget jelent a Monte Carlo integrálás során,
ha az adott eloszlásfüggvény túlzottan összetett, azaz nehezen vehetünk mintákat. Az MCMC módszerek gyökere tulajdonképpen ennek a problémának az áthidalása.
8
A feladatunk tehát egy megfelel®
P (x, y)
átmenetmag találása.
Keressük ezt a
következ® alakban:
P (x, dy) = p(x, y) + r(x)δx ahol
p(x, x) = 0, δx
a Dirac-féle delta függvény, valamint
r(x) = 1 −
R
p(x, y)dy
annak
a valószín¶sége hogy egy lépés során a lánc állapota nem változik. Ez a valószín¶ség nem lehet
0,
hiszen a kifejezésben szerepl® integrál értéke nem feltétlenül
1.
Ha a
p(x, y) függvényre teljesül a megfordíthatóság feltétel, akkor P (x, ·) a π(·) kívánt eloszláshoz fog tartani, ha az iterációk számával tartunk a végtelenbe. Ennek belátásához értékeljük ki az els® egyenl®ség jobb oldalát, el®ször helyettesítsünk be
Z Z
Z P (x, A)π(x)dx =
P -t:
Z p(x, y)dy π(x)dx + r(x)δx (A)π(x)dx =
A cseréljük meg az integrálokat, valamint mivel
Z Z =
δx A
felett pontosan
1,
Z
p(x, y)π(x)dy dy + A
r(x)π(x)dx = A
a megfordíthatósági feltétel miatt
Z Z =
r(x)π(x)dx =
p(y, x)π(y)dx dy + A
r(y)
Z A
de níciója alapján
(1 − r(y))π(y)dy +
=
Z
Z
Z A
r(x)π(x)dx = A
π(y)dy. A
A Metropolis-Hastings algoritmus célja, hogy találjon egy teljesül a fenti tulajdonság.
9
p(x, y) függvényt, melyre
3. fejezet Metropolis-Hastings módszer
3.1.
A Metropolis-Hastings algoritmus
A célunk az lesz, hogy mintákat generáljunk egy
f
egy s¶r¶ségfüggvény,
K
π=
f alakú s¶r¶ségfüggvényb®l, ahol K
pedig a normalizáló konstans.
Tegyük fel, hogy van egy s¶r¶ségfüggvényünk, mellyel generálhatunk jelölteket. Legyen ez
q(x, y),
ahol
R
q(x, y)dy = 1.
x
állapotban vagyunk akkor egy
y
esetén
q(x, y)
y
Ez szemléletesen azt jelenti, hogy ha épp az
q(x, y)-b®l.
értéket generálunk
Ha minden
x
és
kielégítené a megfordíthatósági feltételt, akkor nem is kellene tovább
keresnünk, de sajnos könnyen lehet hogy például
π(x)q(x, y) > π(y)q(y, x). Ez azt eredményezi, hogy a folyamat túl gyakran lép az pedig túl ritkán. Az ötlet ennek kiküszöbölésére, hogy az egy
α(x, y) < 1
valószín¶séggel határozzuk meg. Ez az
x
hogy a folyamat bizonyos esetekben nem lép Így az
x-b®l y -ba (y 6= x)
x
állapotból
y -ba,
fordítva
x állapotból y állapotba lépést
α
sok esetben azt eredményezi,
állapotból
y -ba,
hanem marad
x-ben.
lépés valószín¶sége:
pM H (x, y) ≡ q(x, y)α(x, y). Az egyenl®tlenségb®l adódóan
α(y, x)-et
mivel ez egy valószín¶ség, a fels® korlátja
1.
10
minél nagyobbnak célszer¶ választani és
Viszont
α(x, y)-t
már nem választahatjuk szabadon, mert
pM H (x, y)-nak
ki kell
elégítenie a megfordíthatósági feltételt, így
π(x)q(x, y)α(x, y) = π(y)q(y, x)α(y, x). Ha
α(y, x) = 1,
akkor kifejezés után:
irányba nézne, akkor természetesen
α(x, y) =
π(y)q(y,x) . Ha a relációs jel a másik π(x)q(x,y)
α(x, y)-t választanánk 1-nek és α(y, x)-et fejeznénk
ki. Ezen választás esetén kiegyensúlyoztuk a két oldalt, azaz fennáll az egyenl®ség. Ezek alapján legyen
min π(y)q(y,x) , 1 , π(x)q(x,y) α(x, y) = 1,
ha
π(x)q(x, y) > 0
egyébként.
Az átmenetmag de níciójának teljességéhez szükségünk van még arra a (nem nulla) valószín¶ségre, hogy a folyamat helyben marad, legyen ez:
Z r(x) = 1 −
q(x, y)α(x, y)dy.
Így az Metropolis-Hastings lánc átmenetmagja:
Z PM H (x, y) = q(x, y)α(x, y) + 1 − q(x, y)α(x, y)dy δx . Mivel
pM H (x, y)-t
úgy választottuk, hogy teljesítse a megfordíthatósági feltételt, a
módszer átmenetmagja tartani fog a
π(x)
stacionárius eloszláshoz.
Megjegyzések: (1) Ha egy jelölt elutasításra kerül, akkor az aktuális állapot lesz a soron következ® is egyben. (2) Szerencsére
α(x, y)
kiszámításához nem szükséges a normalizáló konstans isme-
rete, ugyanis megjelenik a számlálóban és a nevez®ben is, így egyszer¶síthetünk vele. (3) Ha szimmetrikus a jelölt generálás, azaz valószín¶sége
π(y)/π(x),
így ha
kerülünk, egyébként pedig
q(x, y) = q(y, x),
π(y) ≥ π(x),
π(y)/π(x)
akkor
1
akkor az állapotváltás
valószín¶séggel
valószín¶séggel.
y
állapotba
Ez szemléletesen azt je-
lenti, hogy ha a generált jelölt értéke nagyobb mint az éppen aktuális, akkor biztos hogy lecseréljük, ha kisebb, akkor pedig
11
π(y)/π(x)
valószín¶séggel.
(4) Az algoritmus alapvet® pontja
q(x, y) megválasztása.
Erre nincs általános eljárás,
de sok praktikus megoldást ismerünk.
Az eljárás összegzése algoritmikus formában:
Válasszunk tetsz®legesen egy FOR
x(0) > 0
kezd®értéket
j = 1, 2, . . . , N
Generáljunk IF
értékeket
u ≤ α(x(j) , y) →
ELSE SET RETURN
y
q(xj , ·)-b®l
SET
és
u
értékeket
U (0, 1)-b®l
x(j+1) = y
x(j+1) = x(j)
{x(1) , x(2) , . . . , x(N ) }
Mint általában az MCMC algoritmusoknál, az értékek a vizsgált
π(x)
eloszlásból
vett mintának tekinthet®ek, ha már a kezd® értéknek semmilyen hatása nincs a folyamatra. Természetesen ehhez szükség van az els® fejezetben részletezett feltételekre, azaz a Markov lánc irreducibilis és aperiodikus mivoltára, ezek viszont nem határozzák meg a konvergencia sebességét.
3.2.
Megvalósítás
Mint a megjegyzéseknél említettük, a megvalósítás szempontjából alapvet® kérdés
q(x, y) megválasztása, melyre sajnos nincs általános eljárás.
A továbbiakban pár gyak-
ran alkalmazott lehet®séget részletezünk.
(1) A jelölt generálás egy lehet®sége ha s¶r¶ségfüggvény. ahol
z
egy
q1
Ezért az
y
q(x, y) = q1 (y − x),
jelöltet az
y = x+z
szerinti valószín¶ségi változó.
ahol
q1 (·)
többváltozós
folyamat szerint választjuk,
Így a jelölt egy jelenlegi és egy
zaj értékének összege, ezt az esetet véletlen bolyongásnak is szokás nevezni. A véletlen mennyiség lehetséges megválasztása például a többváltozós normális vagy t-eloszlás.
12
Amikor
q1
szimmetrikus és
q1 (z) = q1 (−z), α(x, y) = min
akkor az állapotváltás valószín¶sége
π(y) ,1 . π(x)
További fontos kérdés a jelöltgenerálás szórása, mely er®sen befolyásolhatja az algoritmus hatékonyságát, ugyanis hatással van arra, hogy a folyamat mennyire gyakran fogadja el a generált jelöltet. Tegyük fel, hogy a mintákat nagyrészt a módusz körül választjuk és a Markov lánc konvergál. Ha a szórás túl nagy, akkor egyes jelöltek túlságosan távol esnek az aktuális értékt®l és így kis valószín¶séggel cserélünk. A szórás csökkentésével kiküszöbölhetjük ezt a problémát, de túl kicsire sem választhatjuk azt, mert akkor az elfogadott állapotok sorozatának autokorrelációja túl nagy lesz, és így lassú a konvergencia. 2001-ben Roberts és Rosenthal megmutatta, hogy az optimális érték
α = 0.23
körül van, ha a közelíteni kívánt eloszlás azonos peremeloszlások
szorzata. (2) Választhatjuk
q(x, y)-t q2 (y)-nak
is, így az el®z® esettel ellentétben az aktuális
állapottól teljesen függetlenül állítunk jelöltet, viszont ahhoz hasonlóan
q2
lehet
többváltozós normális illetve t-eloszlású, de ebben az esetben meg kell adnunk a várható értéket és a szórást. Ekkor a Metropolis-Hastings hányados
α(x, y) = min
π(y)q2 (x) ,1 . π(x)q2 (y)
Ennek a választásnak egy érdekes esete, ha értéke pontosan
1,
ezáltal a kívánt eloszlásból vesszük a mintákat.
(3) Kézenfekv® megoldás az ismert
π(·)
van rá. Tegyük fel például, hogy vehetünk, pán
v(t)
v(t)-re
q2 (y) = π(y), hiszen ekkor a hányados
eloszlás felhasználása, amikor lehet®ségünk
π(t) ∼ v(t)h(t),
ahol
h(t)
pedig egyenletesen korlátos, ekkor legyen
eloszlásból mintákat
q(x, y) = h(y).
Így csu-
van szükségünk az átmenetvalószín¶ségek meghatározásánal. Ekkor
α(x, y) = min
13
v(y) ,1 . v(x)
4. fejezet Gibbs Sampling
4.1.
A Gibbs mintavételez® eljárás
A Gibbs mintavételez® eljárás a Metropolis-Hastings algoritmus egy speciális esete, így a feladat ugyanaz: hogyan konstruáljunk olyan Markov láncot, amely konvergál a kívánt eloszláshoz.
Segítségével véletlen mintákat generálhatunk egy többváltozós
valószín¶ségeloszlásból anélkül, hogy ismernénk a s¶r¶ségfüggvényt. Ezeket a mintákat a kés®bbiekben használhatjuk az együttes eloszlás, egy változó peremeloszlásának vagy a változók egy részhalmazának együttes eloszlásának közelítéséhez. Tegyük fel hogy adott egy
f (x) = (x, y1 , . . . , yp )
Z
együttes s¶r¶ségfüggvény:
Z
f (x) =
...
f (x, y1 , y2 , . . . , yp )dy1 dy2 . . . dyp .
Számos esetben nagyon nehéz egyszer¶en kiszámolni juk a Gibbs mintavételez®t, hogy generáljuk egy
f (x)-et,
ilyenkor használhat-
X1 , . . . , Xm ∼ f (x)
sorozatot
f (x)
ismerete nélkül. Ha elég nagy mintát generálunk, a várható érték, a szórás, vagy még akár a s¶r¶ségfüggvény isg meghatározható egy adott hibahatáron belül. Az eljárás bemutatásához el®ször tekintsünk egy kétváltozós esetet. valószín¶ségi változó párok
(X, Y ).
generálna mintát, az (ismert)
A Gibbs mintavételez® ahelyett, hogy
f (x | y)
f (y | x)
és
0
0
0
0
0
0
0
Y0 , X0 , Y1 , X1 , Y2 , X2 , . . . , Yk , Xk .
14
f (x)-b®l
feltételes eloszlásokat használja,
mégpedig valószín¶ségi változók egy Gibbs sorozatát gyártja le:
0
Legyenek a
Az
0
0
Y0 = y 0
kezd®érték ismeretében a további értékeket iteratívan határozzuk meg
a következ® generálási szabály szerint:
0
0
0
Xj ∼ f (x | Yj = yj ) 0
0
0
Yj+1 ∼ f (y | Xj = xj ). Az így generált sorozatot Gibbs mintának nevezzük. Megfel®len általános feltételek mellett
0
Xk
Gibbs minta legutolsó eleme
4.2.
f (x)-hez,
eloszlása konvergálni fog
0
0
X k = xk ,
ha
k → ∞.
Így ha
k
elég nagy, akkor a
f (x)-b®l.
egy meg gyelt érték
Konvergencia
Egyáltalán nem nyilvánvaló, hogy egy
f (x)
eloszlású valószín¶ségi változó el®állítható
a Gibbs minta segítségével, illetve hogy ez konvergál. Ez az iteráció Markov tulajdonságán múlik. El®ször nézzük a következ® esetet: legyenek
X
és
Y
Bernoulli peremel-
oszlású valószín¶ségi változók a következ® együttes eloszlással:
P (X = 0, Y = 0) = p1 ,
P (X = 1, Y = 0) = p2 ,
P (X = 0, Y = 1) = p3 ,
P (X = 1, Y = 1) = p4 ,
ahol
pi ≥ 0,
és
P
pi = 1.
i Mátrixos alakban az együttes eloszlással:
fx,y (0, 0) fx,y (1, 0) fx,y (0, 1) fx,y (1, 1) Ez alapján
x
! =
p1 p2 p3 p4
! .
peremeloszlása:
fx = (fx (0), fx (1)) = [(p1 + p3 ), (p2 + p4 )]. Ezek alapján a feltételes eloszlásokat könnyen számolhatjuk, mátrix alakban:
p3 p1 +p3 p4 p2 +p4
!
Ay|x =
p1 p1 +p3 p2 p2 +p4
p2 p1 +p2 p4 p3 +p4
!
Ax|y =
p1 p1 +p2 p3 p3 +p4
,
valamint
15
,
ahol
Ay|x
az
Y -nak (X = x)-re
vonatkozó feltételes valószín¶sége és fordítva.
Ha ezekre vonatkozóan el®állítjuk a Gibbs mintát, akkor egy kapni.
Legyen
lépés és
Ay|x
Ax|y
az az átmenetmátrix, amely az
y -okból x-ekbe
pedig az
Tekintsük csupán
X
x
0 − 1 sorozatot fogunk y
állapotokból az
állapotokba
lépés átmenetvalószín¶ségeit tartalmazza.
peremeloszlásának közelítését, ekkor a Gibbs minta
van szükségünk. Ahhoz, hogy
0
0
X0 -ból X1 -et 0
meghatározzuk, szükségünk van
0
0
X0 → Y1 → X1
az iterációval kapott sorozatunk
0
Xi
és
0
0
X0 → X1
elemeire 0
Y1 -re,
így
Markov lánc lesz a
következ® átmenetvalószín¶séggel:
0
0
P (X1 = x1 | X0 = x0 ) =
X
0
0
0
0
P (X1 = x1 | Y1 = y) · P (Y1 = y | X0 = x0 ).
y Általánosan, az
X
0
Ax|x
sorozat átmenetmátrixa,
az következ® módon határozható
meg:
Ax|x = Ay|x Ax|y . Ennek segítségével bármely alapján
0
0
P (Xk = xk | X0 = x0 )
k -ra
meghatározhatjuk
pontosan
(Ax|x )
k
0
Xk -t,
. Továbbá, ha
mert az átmenetmátrix 0
Xk
peremeloszlását
fk = [fk (0), fk (1)] k -ra
alakban írjuk fel, akkor akkor minden
fk = f0 (Ax|x )k = (f0 (Ax|x )k−1 )Ax|x = fk−1 Ax|x . Ha
Ax|x
minden eleme pozitív, akkor a fentiek alapján minden
n¶ség esetén, ha
k → ∞ , fk
konvergál az
f
f0
kezdeti valószí-
stacionárius elolszláshoz, melyre teljesül,
hogy
f Ax|x = f. Emiatt
f
az
X
változó peremeloszlása. Ez alapján kell®en sokszor iterálva tetsz®-
legesen közel kerülhetünk a
X
peremeloszlásához, azaz
k→∞
esetén
0
Xk
közelít®leg
fx . Az eddig tárgyalt
2 × 2-es eset könnyen átvihet® bármilyen n × m-es együttes elosz-
lásra. Hasonlóan de niálhatjuk az
n×n-es Ax|x átmenetmátrixot, melynek stacionárius 16
eloszlása
X
peremeloszlása lesz, viszont folytonos
van szükségünk.
0
X
és
Y
esetén kisebb módosításokra
0
X1 X0 -ra vonatkozó feltételes eloszlását a következ® alakban írjuk fel: Z fX 0 |X 0 (x1 | x0 ) = 1
fX 0 |Y 0 (x1 | y)fY 0 |X 0 (y | x0 )dy,
0
valamint a feltételes eloszlások:
1
0
1
1
0
0
0
0
0
0
X1 | X0 , X2 | X3 , X1 | X0 ,
....
Ax|x -hez
hasonlóan
az átmenetmátrixnak ki kell elégítenie a következ® egyenl®séget:
Z f Itt
fX 0 |X 0 k
0 0 Xk |X0
(x | x0 ) =
fX 0 |X 0 (x | t)fX 0 k
k−1
0 k−1 |X0
(t | x0 )dt.
személeletesen az átmenet eloszlása. Hasonlóan a Bernoulli eloszlású
k−1
változókhoz, ha
k → ∞ , fk
szintén a stacionárius eloszláshoz tart, ami pedig
X
pe-
remeloszlása lesz.
4.3.
Két- és többváltozós eset
El®ször tegyük fel, hogy adottak és
fY |X (y | x)
X
Y
és
valószín¶ségi változók, valamint
feltételes eloszlásfüggvények. Határozzuk meg az
fX (x)
fX|Y (x | y)
peremeloszlást.
De níció szerint
Z fX (x) = ahol
fXY (x, y)
fXY (x, y)dy,
az ismeretlen együttes eloszlás. Mivel
fXY (x, y) = fX|Y (x | y)fY (y),
behelyettesítve
Z fX (x) = Ha
fY (y)-t
fX|Y (x | y)fY (y)dy.
hasonlóan helyettesítjük
Z fX (x) =
Z fX|Y (x | y)
fY |X (y | t)fX (t)dtdy =
Z Z fX|Y (x | y)fY |X (y | t)dy fX (t)dt = h(x, t)fX (t)dt,
ahol
h(x, t) = fX|Y (x | y)fY |X (y | t)dy .
egyértelm¶ megoldása.
17
Ennek az integrálegyenletnek
fX (x)
Ezáltal ha
k → ∞, fX 0 |X 0 (x | x0 ) → fX (x), 0
k
valamint
fX 0 |X 0 (x | t) → h(x, t). k
Sajnos annak ellenére, hogy
k−1
X
és
Y
együttes eloszlása meghatározza a feltételes-
és peremeloszlásokat, a megfelel® feltételes eloszlások nem határozzák meg a megfelel® peremeloszlásokat és ezáltal az együttes eloszlást sem. Ha több változónk van, a helyzet bonyolultabbá válik, hiszen nem minden feltételesés peremeloszlás pár határozza meg az együttes eloszlást, ezáltal többféle módon is felírhatjuk a fenti integrálegyenletet, így a feltételes eloszlások különböz® halmazait használhatjuk, hogy meghatározzuk a kívánt peremeloszlást. Tegyük fel, hogy
fX (x)-et
keressük, valamint
A fenti integrálegyenlethez hasonlóan, ha az
X ,Y
(Y, Z)
és
Z
a valószín¶ségi változók.
párt egy valószín¶ségi változóként
értelmezzük:
Z Z
fX|Y Z (x | y, z)fY Z|X (y, z | t)dydz fX (t)dt.
fX (x) = fX|Y Z
konvergál.
fY |XZ
és
fY Z|X
és
váltakozása szintén el®állít egy sorozatot, mely az
fX (x) eloszláshoz
Ezzel szemben a Gibbs mintavételez® iteratívan venne mintákat
fZ|XY -ból,
így a
k -adik
fX|Y Z ,
iteráció után:
0
0
0
0
0
Xj ∼ f (x | Yj = yj , Zj = zj ) 0
0
0
0
0
Yj+1 ∼ f (y | Xj = xj , Zj = zj ) 0
0
0
0
0
Zj+1 ∼ f (z | Xj = xj , Yj+1 = yj+1 ). Ez az iteráció a következ® sorozatot gyártja le:
0
0
0
0
0
0
0
0
0
Y0 , Z0 , X0 , Y1 Z1 , X1 , Y2 , Z2 , X2 , . . . és ha
k -t
megfelel®en nagynak választjuk,
értéke.
18
0
0
X k = xk
közelít®leg
f (x)
egy meg gyelt
Mivel a Gibbs mintavételez® felhasználja az összes egyváltozós eloszlást, ezek meghatározzák az együttes és az összes feltételes eloszlást, ezért az integrálegyenlet megoldható az iterációs séma segítségével.
4.4.
Megvalósítás
A megvalósítás során az alapvet® kérdés az, hogy milyen nagynak válasszuk a
k -t, hogy
az eltérés egy bizonyos határon belül maradjon. A legáltalánosabb módszer a Gibbs minta konvergenciájának vizsgálata bizonyos tényez®kt®l függ®en. Generálhatunk darab független Gibbs mintát és választhatjuk úgy a
m
k -t, hogy ezek a minták valamilyen
statisztikai próba szerint azonos eloszlásúak legyenek, vagy vizsgálhatjuk a közelít® és a közelített eloszlások közötti különbségek változását. A megfelel® is számos kutatás tárgyát képezi.
19
k
megválasztása ma
5. fejezet Slice sampling
A feladatunk hasonló az eddigiekhez, szeretnénk egy mintát vételezni, ahol
x ∈ Rn .
p(x) =
A slice sampling ötlete, hogy ezt megtehetjük úgy,
hogy egyenletes eloszlású mintákat veszünk az
f (x)
gra konja alatti
tartományból. Az ötlet formalizálásához vezessünk be egy
együttes eloszlásfüggvényét úgy, hogy egyenletes legyen az
f (x)}
által meghatározott tartományon. Ezáltal
( p(x, y) = Z=
R
f (x)dx.
Az
x
(x, y)
1/Z
ha
0
egyébként,
(n + 1)
dimenziós
y segédváltozót és de niáljuk
x és y
ahol
f (x) alakú eloszlásból K
U = {(x, y) : 0 < y <
együttes eloszlása:
0 < y < f (x)
peremeloszlása ekkor
Z p(x) =
f (x)
1/Zdy = f (x)/Z. 0
Így
(x, y)-ra
hagyjuk
y -t.
vonatkozó mintákat fogunk gyártani, majd egyszer¶en gyelmen kívül
Vegyük észre, hogy ezáltal a slice sampling egy speciális esete a Gibbs
mintavételez®nek, ugyanis ha
(X, Y )
az együttes eloszlás, akkor az
Y | X = x ∼ U (0, f (x)), X | y ∼ U (S),
ahol
S = {x : y < f (x)}
egy Gibbs sorozat.
20
5.1.
Egyváltozós eset
A legegyszer¶bb eset az, ha a keresett eloszlásfüggvény egyváltozós. A továbbiakban többváltozós eloszlásfüggvényeket fogunk tekinteni, de az
n-dimenziós mintát koordiná-
xi
koordináta esetében meg kell
tánkénti mintavételezéssel fogjuk legyártani. Minden határoznunk
fi (xi ) értéket, mely a K
Az egyszer¶ség kedvéért
n-dimenziós
x
normalizáló konstanstól eltekintve
jelöljön egy
n-dimenziós
eloszlásfüggvényt. A következ® eljárás az aktuális
1 Generáljunk egy
y ∼ U (0, f (x0 ))
tározunk egy vízszintes
S
f (x)
vektort,
x0
p(xi | {xj }j6=i ). pedig jelöle az
vektort
x1 -re
cseréli:
egyenletes eloszlású valós számot, ezzel megha-
szeletet:
S = {x : y < f (x)}. x0
mindig
S -en
belül
lesz. 2 Keressünk egy
E = (L, R)
intervallumt
x0
körül.
Ha
x0 ∈ S ,
akkor
x0
egy
környezete mindig tartalmazni fogja a szelet valamekkora részét. 3 Generáljunk egy
x1
egyenletes eloszlású véletlen pontot az
Az els® lépésben tehát generálunk egy szeletet.
y
E
intervallumon.
értéket és meghatározzuk a vízszintes
A második és harmadik lépést különböz® módokon is megvalósíthatjuk, de
ezeknek természetesen egy megfelel® Markov láncot kell eredményezniük. Lehet®ségek alkalmas
E
intervallum keresésére:
1 A lehet® legjobb megoldás az lenne, ha
L = inf(S)
és
R = sup(S)
lenne, de ez
általában nem kivitelezhet®. 2 Ha a lehetséges
x
értékek korlátosak, akkor
E -t
ehhez az korláthoz hasonlóan
állítjuk be. 3 Ha van egy
w
becslésünk az
hosszúságú intervallumot
x0
S
méretére, véletlenszer¶en kiválasztunk egy
körül, majd addig növeljünk, amíg le nem fedi
w
S -t,
például a stepping out eljárással. 4 Hasonlóan, véletlenszer¶en kiválasztunk egy
w hosszúságú intervallumot x0 körül,
majd növeljük a duplázó eljárással.
Mindegyik eljárás esetén meg kell határoznunk egy jelölteket tartalmazza, a következ® szerint:
21
A halmazt, mely az elfogadható
A = {x : x ∈ S ∩ E ∧ P (E | x) = P (E | x0 )}. Ez természetes elvárás, ugyanis az elfogadható jelöltnek a szelet és a választott intervallum metszetében kell lennie, illetve az aktuális és a soron következ® állapotot
E -ben
feltéve ugyanolyan valószín¶séggel kell teljesülnie, hogy Könnyen láthatjuk, hogy az
1.
és
2.
esetben
A = S.
Az
vagyunk.
1.
esetben
E -t
pont így
állítottuk be, de ez csak akkor hatékony, ha valamilyen módszer segítségével el® tudjuk állítani
y = f (x)
lehetséges
x
összes megoldását.
A
2.
lehet®séget könny¶ megvalósítani, ha a
jelöltek halmaza korlátos, viszont nem túl hatékony, ha a szelet általában
sokkal kisebb ennél a korlátnál. A
3.
pontban említett stepping out eljárás algoritmusa:
f
Input: natkozó
w
függvény,
x0
becslés és az
Output:
(L, R)
Generáljunk
y
pont, a szeletet meghatározó
m
érték, a szelet maximális
mw
érték, a szelet méretére vo-
méretére vonatkozólag.
intervallum.
U ∼ U (0, 1)
és
V ∼ U (0, 1)
egyenletes eloszlású véletlen számokat
L ← x0 − w · U ; R ← L + w; J ← [m · U ]; K ← (m − 1) − J REPEAT WHILE
(J > 0
AND
y < f (L)):
L ← L − w; J ← J − 1 REPEAT WHILE
(K > 0
AND
y < f (R)):
R ← R + w; K ← K − 1
5.1. Megjegyzés. Itt [m · V ] az m · V egész részét jelöli. A stepping out eljárást abban az esetben alkalmazhatjuk, ha tudunk mondani egy becslést az
S
szelet jellemz® méretére. Ekkor
fels® korlátja, ahol
w
+
m∈Z
mw
. Speciális esetben
w
a keresett intervallum hosszának egy
m = 1,
ekkor az intervallumot mindig
hosszúságúra állítjuk be, így nem szükséges meghatározni
f
értékét az intervallum
végpontjaiban. Nagyon fontos, hogy az intervallum kezdeti elhelyezkedése, valamint az intervallum bal és jobb irányba történ® növelésének elosztása véletlenszer¶ legyen, mert ez biztosítja, hogy a végeredményben kapott intervallum bármelyik kezdeti pontból el®állítható legyen.
22
x0 ∈ S ∩ E
A
4.
pontban említett doubling eljárás algoritmusa:
f
Input: natkozó
w
függvény,
becslés és a
(L, R)
Output:
x0 p
pont, a szeletet meghatározó
érték, a szelet maximális
2p w
y
érték, a szelet méretére vo-
méretére vonatkozólag.
intervallum.
Generáljunk egy
U ∼ U (0, 1)
egyenletes eloszlású véletlenszámot.
L ← x0 − w · U ; R ← L + w ; K ← p REPEAT WHILE
(K > 0
Generáljunk egy IF
V <
{y < f (L)
V ∼ U (0, 1)
OR
y < f (R)}:
egyenletes eloszlású véletlenszámot.
1 2
THEN ELSE
AND
L ← L − (R − L)
R ← R + (R − L)
K ←K −1 Szemléletesen minden lépésben egy intervallumot gyártunk le, melynek hossza az el®z® intervallum hosszának kétszerese. Ezt addig ismételjük, amíg a soron következ® intervallum mindkét végpontja a szeleten kívül esik, vagy az el®re meghatározott lépésszámot el nem érjük.
Fontos megjegyezni, hogy minden lépésben véletlenszer¶en
döntjük el, hogy melyik irányban növeljük az intervallumot (emiatt az egyik végpont változatlan marad), attól függetlenül, hogy az adott irány szerinti végpont a szeleten belül vagy azon kívül helyezkedik el. A véletlenszer¶ség a doubling eljáráshoz hasonlóan alapvet®, mert itt is ez garantálja, hogy az intervallum más kezdeti
x0 ∈ S ∩ E
is el®állítható legyen. Sajnos ebben az esetben viszont nem minden
pontból
x∈S∩E
pontra
teljesül, hogy onnan indulva el®állíthatjuk ugyanazt az intervallumot, ezért szükséges ellen®riznünk, hogy a következ® állapot
A-beli
vagy sem.
A stepping out eljárással szemben ez a módszer gyorsan növeli a kezdeti intervallum méretét, ami nagyon hatékony, ha a szelet becsült
w
értéke kicsi, ugyanis hasonlóan
kicsi kezdeti intervallummal indulunk és mivel exponenciálisan n®, hamar lefedjük a szeletet.
23
5.1.1. ábra: A stepping out, shrinkage és doubling eljárások. Az (a) részben meghatározzuk körül, majd
y -t. A (b) részben egy w hosszúságú intervallumot véletlenszer¶en elhelyezünk x0
w méret¶ intervallumok hozzáadásával addig növeljük, amíg az intervallum mindkét vége a szeleten kívül
helyezkedik el. A (c) részben addig csökkentjük az intervallum méretét, amíg az új
x1 pont az szeleten belül nem
helyezkedik el. A (d) részben kétszer duplázzuk az intervallum méretét, ekkor már mindkét vége a szeleten kívül esik.
Ha megtaláltuk a megfelel® intervallumot, akkor még generálnunk kell egy véletlen számot, mely a szeletre esik:
E
intervallum
E
kiinduló in-
tervallumon, majd a shrinkage módszerrel folyamatosan csökkentjük
E -t egészen
1 Egymás után generáljunk egyenletes eloszlású véletlen számokat az egész addig, amíg nem kapunk egy
A-belit.
2 Egymás után generáljunk egyenletes eloszlású véletlen számokat az
addig, amíg nem kapunk egy
A-belit.
24
A
2.
pontban említett shrinkage eljárás algoritmusa:
Input:
f
függvény,
x1
Output:
x0
pont, a szeletet meghatározó
y
érték,
(L, R)
intervallum.
pont.
0
0
L ← L; R ← R REPEAT
U ∼ U (0, 1)
Generáljunk egy 0
0
egyenletes eloszlású véletlen számot.
0
x1 ← L + U · (R − L ) IF
y < f (x1 )
AND Accept(x1 )
THEN EXIT LOOP IF
x1 < x0
THEN ELSE
0
L ← x1 0
R ← x1
5.2. Megjegyzés. Itt az Accept(x1 ) egy tesztet jelöl, arra vonatkozólag hogy x1 ∈ S∩E
vagy sem.
Az algoritmus futása során minden egyes lépésben, amikor olyan pontot kapunk, mely nincs benne az
S ∩E
halmazban, az aktuális intervallum összezsugorodik, mégpe-
dig úgy, hogy az aktuális pont lesz az egyik végpont. Mivel intervallumban lesz
A-beli
x0
mindig
A-beli,
a kapott
pont, így a folyamat biztosan végetért. Az intervallum fo-
lyamatos csökkentése miatt a sorsolt pontok várható száma kicsi, ezért ez a módszer általában jól használható. Az Accept(x1 ) teszt: Input:
f
függvény,
tére vonatkozó Output: 0
x1
x0
pont,
x1
jelölt, a szeletet meghatározó
érték, a szelet mére-
w becslés, w becsléssel és a doubling eljárással talált (L, R) intervallum. elfogadható következ® pontnak vagy sem. 0
L ← L; R ← R D←
y
FALSE
25
REPEAT WHILE 0
0
0
R − L > 1, 1 · w
0
M ← (L + R )/2 IF
{x0 < M
THEN IF
D←
x1 ≥ M }
OR
{x0 ≥ M
AND
x1 < M }
TRUE
x1 < M
THEN ELSE IF
AND
D
0
R ←M 0
L ←M
AND
0
y ≥ f (L )
AND
0
y ≥ f (R )
THEN nem fogadjuk el az új pontot Az új pontot akkor fogadjuk el, ha a fenti eljárás során nem lett elutasítva
5.3. Megjegyzés. Az 1, 1 szorzó a numerikus számítás során létrejöv® kerekítési hiba
kiküszöbölésére szolgál.
Az eljárás ellen®rzi a doubling módszerrel legyártott intervallumokat, mégpedig vizsgálja hogy azok végei a szeleten kívül esnek-e vagy sem, ugyanis el®bbi az algoritmus korábbi befejezését jelentené. Ha a megfelel® intervallumot az akkor
5.2.
A = S ∩ E.
A
(4)
(1) − (2) − (3)
eljárással lehet, hogy az
módszerek valamelyikével kerestük,
A
az
S∩E
egy részhalmaza.
Többváltozós eset
Ahelyett, a többváltozós eloszlás minden
xi
koordinátájára alkalmaznánk az eddig is-
mertett módszerek egyikét, a megfelel® fogalmak
n-dimenziós változatának használatá-
val értelmezhetjük a módszereket magára az
n-dimenziós eloszlásra.
E = (L, R)
n-dimenziós
intervallum helyett vegyünk egy
téglát:
H = {x : Li < xi < Ri ∀ i = 1, . . . , n}.
26
Az eddig használt
Az
x1 = (x1,1 , x1,2 , . . . , x1,n ) következ® állapot megtalálása az aktuális x0 = (x0,1 , x0,2 , . . . , x0,n )
mellett hasonló az egyváltozós esethez:
1 Generáljunk egy határozunk egy 2 Keressünk egy
y ∼ U (0, f (x0 ))
S
egyenletes eloszlású valós számot, ezáltal meg-
S = {x : y < f (x)}.
szeletet:
H = (L1 , R1 ) × (L2 , R2 ) × . . . × (Ln , Rn ) n-dimenziós
téglát
x0
körül. 3 Generáljunk egy egyenletes eloszlású
x1 -et az S
szelet azon részén, ami a
H
téglán
belül van.
Megfelel®
H
tégla kereséséhez az egyváltozós esethez hasonlóan ideális lenne az
tartalmazó legkisebb korlátos, akkor
H -t
H
megtalálása, de ez nehezen kivitelezhet®. Ha minden változó
a korlátokhoz hasonlóan állíthatjuk be, de ez mint az egyváltozós
esetben, itt sem hatékony, ha az egy
wi
becslésünk az
w = [wi ]
téglát az
x0
S
S
szelet sokkal kisebb. Ha minden tengely mentén van
lehetséges méretére, akkor véletlenszer¶en választhatunk egy
pont körül, majd növelhetjük valamilyen eljárással. A shrinkage
eljárás során felmerül az a nehézség, hogy mivel az nagy
n esetén túl sok tesztet kéne végeznünk.
igényes lehet az
n
S -et
n-dimenziós téglának 2n
csúcsa van,
A stepping out módszer is túlzottan id®-
lehetséges irány miatt. A doubling eljárás megvalósításánál szintén
akkor állunk meg, amikor egy egyenletes eloszlás szerinti véletlenszám az
S
szeleten
kívül esik.
5.3.
Egyenletes eloszlás egy kétdimenziós tartományon
Mint az eddigi módszereknél, a slice sampling alkalmazhatósága azon múlik, hogy a generált Markov lánc stacionárius eloszlása megegyezik-e a szóban forgó eloszlásfüggvénnyel.
Ezt nem láttuk be a fenti esetekben, viszont adunk egy szemléletes és
heurisztikus bizonyítást a kétdimenziós esetre. Legyen
T
egy tartomány a síkon, melynek határa nullmérték¶. Generáljunk ezen a
tartományon egyenletes eloszlású
v = (x, y)
koordinátájú véletlen pontot a következ®
átmenetszabály szerint:
1 Az els® lépésben válasszunk tetsz®legesen egy vízszintes egyenest, melynek van olyan szakasza, ami a tartományon belül van. Állítsuk be az
27
y
koordinátát
y=
y1 -re.
Válasszuk ki tetsz®legesen az egyik
végpontja legyen
L1 ,
a jobb pedig
R1 .
eloszlású véletlenszámot. Állítsuk be 2 Az
(x1 , y1 )
mazó,
m1
l1
hosszúságú szakaszt, melynek bal
Generáljunk egy
u1 = U (0, l1 )
x1 = L1 + u1 -t.
ponton átmen® függ®leges egyenes mentén legyen az
D1 ,
hosszúságú szakasz alsó végpontja legyen
neráljunk egy
u2 = U (0, m1 )
egyenletes
(x1 , y1 )-t
a fels® pedig
egyenletes eloszlású véletlenszámot.
tartal-
U1 .
Ge-
Állítsuk be
y2 = D1 + u2 -t. 3 Az
(x1 , y2 )
mazó,
m2
ponton átmen® vízszintes egyenes mentén legyen az
L2 ,
hosszúságú szakasz bal végpontja legyen
neráljunk egy
u3 = U (0, m2 )
(x1 , y2 )-t
a jobb pedig
egyenletes eloszlású véletlenszámot.
tartal-
R2 .
Ge-
Állítsuk be
x2 = L2 + u3 -t. 4 A (2)-(3) lépéseket ismételjük egymás után.
Tehát minden páratlan lépésben vízszintesen választunk egy új pontot, minden párosban pedig függ®legesen.
Így egy Markov láncot kapunk, melyre teljesülnek a
megfelel® feltételek és stacionárius eloszlása egyenletes.
Lássuk be ezt az állítást a
Markov láncokra vonatkozó feltételek felhasználása nélkül. Fedjük le a tartományt egy
q lépéstávolságú rácshálóval és tekintsük a rácsháló által
meghatározott négyzeteket, ezek közül is azokat, melyekre van olyan az
(x, y)
(x, y) ∈ T ,
hogy
pont a négyzet bels® pontja. Ezen négyzetek középpontjait kössük össze az
oldallapjaikkal szomszédos négyzetek középpontjaival. Bolyongjunk (a paramétereket diszkrét esetben értelmezve) a négyzeteken a fenti átmenetszabály szerint.
Világos,
hogy ha csak két szomszédos négyzetünk van (akár vízszintesen, akár függ®legesen), mindkett®ben
1 1 - valószín¶séggel vagyunk, a teljes rácsháló pedig ilyen egymás melletti 2 2
négyzetek halmaza. Azon pontok esetében, melyeket két négyzet csúcsa fed csupán, egészítsük ki a lefedést két tovább négyzettel, így az ilyen pontoknak lesz olyan sugarú környezete, melyben minden pont valamely négyzet bels® pontja. Így ennek a láncnak a stacionárius eloszlása egyenletes lesz és ha a felosztást minden határon túlmen®en nomítjuk (q
→ 0),
akkor az egész
T
tartományon egyenletes az eloszlás.
A fenti eset az általánosság csekély korlátossága mellett érvényes, de nyilvánvalóan kiterjeszthetjük általánosabb idomokra is.
28
6. fejezet Perfect sampling
6.1.
Az alapötlet
Az eddig tárgyalt módszereknél általános problémát jelentett, hogy milyen nagynak válasszuk
t-t,
hogy
Xt
közelít®leg egy minta legyen a
szempontjából további nehézséget jelent
Xt
és
Xt+1
π
eloszlásból.
A megvalósítás
korrelációjának meghatározása,
mely a közelítés szórásához szükséges. A Perfect sampling mindkét problémát orvosolja azáltal, hogy pontosan a
π
eloszlásból gyárt független mintákat. Az ötlet az MCMC
módszerekre általánosan jellemz® probléma kiküszöbölésén nyugszik: ha megfelel®en hosszú (t
→ ∞)
ideig futtatjuk a láncot, akkor megkapjuk a keresett eloszlást, de
a számítási kapacitás és a tárhely korlátossága miatt meg kell állnunk egy bizonyos id®pillanatban. A Perfect sampler véges sok lépésben legyártja pontosan azt a Markov láncot, melynek végtelen sok lépést kellene tennie azáltal, hogy lecseréljük a kezdeti id®pillanatot eloszlású a
−∞-re, ∞-t
t=0
pedig
0-ra.
Így a futtatási id®
t = ∞,
id®pillanatban és ekkor pontosan a keresett
π
a lánc stacionárius
eloszlásból vételezzük
a mintát. El®ször is egy olyan
Xt -t kell keresnünk, mely független a kiindulási állapottól.
több Markov lánc egyidej¶ futtatásával fogjuk elérni. Tegyük fel, hogy az nek
I
Ezt
állapottér-
k darab állapota van és indítsunk a t = 0 id®pillanatban k láncot minden különböz®
állapotból, ezek lesznek a párhuzamos láncaink. Ezeket kapcsoljuk össze egy netszabállyal és
Ut
véletlen számokkal, ahol
φ
függvényét. Minden lánc esetében ugyanazt a használjuk, ahol általában
Ut+1 ∼ U (0, 1)
meghatározza
mint
Xt
átme-
és
Ut+1
φ függvényt és . . ., Ut , Ut+1 , . . . számokat
valamint
29
Xt+1 -et,
φ
Xt+1 = φ(xt , ut+1 ) = FX−1t+1 |xt (ut+1 )
az átmenetmag és az állapottér lineáris rendezése szerinti függvénye. Jelölje
X s,j ≡ Xts,j t≥s
Xt+1 | xt
azt a Markov láncot, mely a
j
inverz eloszlás-
állapotból indul az
s
id®pillanatban.
6.1. De níció. Azt mondjuk, hogy a Markov láncok egyesültek a t id®pillanatban, ha Xts,1 = Xts,2 = . . . = Xts,k .
6.1. Tétel. Tegyük fel, hogy X s,1 , X s,2 ,. . .,X s,k összekapcsolt Markov láncok, ahol
1.
X s,j
a j -edik állapotból az s id®pillanatban indul, és
2.
s,j Xt+1 = φ(Xts,j , Ut+1 ),
ahol az Ut-k páronként függetlenek.
Ekkor 1. Az egyesülés T id®pillanata egy valószín¶ségi változó, ami csak U1, U2,. . . függvénye, valamint 2. az egyesülés állapota az XT valószín¶ségi változó független a kezdeti állapottól. Az els® állítás azonnal következik a következménye, hogy
XT
állítás jelentése, hogy
T
nem függ
φ
j -t®l
konstrukciójából, a második pedig annak a
U1 ,. . .,UT
mivel csak
az az id®pillanat, amikor
XT
lási állapottól, de sajnos azt nem mondhatjuk, hogy eloszlásból. mivel
6.2.
T
Ha
T
∗
egy x id®pillanat és
egy valószín¶ségi változó,
XT ∗
függvénye. A második
már teljesen független a kiindu-
XT
egy minta a
π
s,j független Xs -t®l, akkor
stacionárius
XT ∗ ∼ π ,
de
XT π .
A Coupling From The Past (CFTP) algoritmus
Propp és Wilson nevéhez köt®dik a CFTP algoritmus, mely x futtatási id®ben kihasználja a láncok egyesülését, ezáltal olyan valószín¶ségi változót gyárt le, melynek pontosan a keresett
π
az eloszlásfüggvénye.
Tegyük fel, hogy egy láncot a állapotból. Ekkor a egy minta a
π
t=0
t = −∞
id®pillanatban indítottunk valamilyen
X−∞
id®pillanatban (a végtelen sok lépés következtében) az
X0
stacionárius eloszlásból. Az implementálás az egyesülés megfontolásán
alapszik: el®ször keressünk egy olyan tehát az egyesülés megtörtént a
−T
−T
és a
-t, amire
0
X0 ≡ X0−T,j
nem függ
−T,j X−T -t®l,
id®pillanat között, majd minden lehetséges
30
állapotból futtatunk egy láncot a
X0 -t.
t = −T
A feladatunk tehát egy megfelel®
id®pillanattól
−T
és
X0
t = 0-ig,
ezáltal megkapjuk
keresése, mely a következ® módon
történik:
1. Indítsuk az
X −1,1 ,X −1,2 ,. . .,X −1,k
állapottér elemszáma 2. Az
k)
és generáljunk egy
−1,j X0−1,j = φ(X−1 , U0 )
a láncok egyesültek a
láncokat a
t = −1
id®pillanatban (ahol az
U0 -t.
átmenetszabály segítségével tegyünk egy lépést.
t = 0
pillanatban, akkor
T = −1
és
X0
Ha
egy minta az
π
eloszlásból. 3. Egyébként indítsuk az
X −2,1 ,X −2,1 ,. . .,X −2,k
−2,j generáljunk egy U−1 -et, majd az X−1
=
láncokat a
t = −2
−2,j φ(X−2 , U−1 ) és
X0−2,j
pillanatban és
−2,j = φ(X−1 , U0 )
átmenetszabályok segítségével tegyünk két lépést. Ha a láncok egyesültek a pillanatban, akkor
T = −2
4. Egyébként vegyük sorban a
és
X0
egy minta a
π
t=0
eloszlásból.
t = −3,−4,. . . pillanatokat és alkalmazzuk az eljárást
iteratívan, amíg a láncok nem egyesülnek a
t = 0-ban.
6.1. Megjegyzés. Fontos, hogy a harmadik lépésben ugyanazt az U0 -at használjuk,
mint az els®ben.
A fenti algoritmus legyárt egy valószín¶ségi változót, melynek eloszlása pontosan megegyezik a Markov lánc eloszlásával. Az eljárással megkeresett
T
és
X0
viszont nem
független változók.
6.3.
Fontainebleau lehulló levelei
Egy látványos geometriai alkalmazása a CFTP módszernek a fákról lehulló levelek elhelyezkedésének modellezése, melyet Fontainebleau fái ihlettek. A modell egy véletlenszer¶en elkészített mozaik kialakulásának folyamatát írja le, melyben egy üveglapra hulló levelek által alkotott képpel közelítjük a stacionárius eloszlást. A levelek lehullásának folyamatát jól modellezi egy Markov lánc, melynek elemei a különböz® minták, melyeket a lehulló levelek alkotnak az üveglapon.
David Wilson nevéhez köt®dik a
CFTP algoritmus egy módosítása, mely lépésenként építi fel a moziakot, visszakövetés nélkül, a végs® képet pedig fokozatosan éri annak a ténynek a felhasználásával,
31
hogy minden újabb lehulló levél a kés®bbiekben csökkenti a kép módosulását az adott területen.
6.3.1.ábra: Fontainebleau lehulló levelei In: [12]
Tekintsünk a levelek alkotta képre az üveglap aljáról, ahelyett hogy felülr®l szemlélnénk. Így abban a pillanatban megsz¶nik változni a kép, amint az üveglap minden egyes pontja le van fedve egy levél által, ezért ekkor az aktuális állapot egy minta a keresett stacionárius eloszlásból.
6.4.
Sandwiching
A CFTP algoritmus kiküszöböli a futási id® hosszának problémáját, de a gyakorlatban nehezen alkalmazható, ha az
I
állapottér
k
elemszáma túl nagy. A sandwiching
ötlete a leginkább elterjedt és talán a leghatékonyabb eljárás a CFTP gyakorlati alkalmazásához nagy állapottér esetén. Ezt azonban csak olyan Markov láncok esetén haszhálhatjuk, melyek kielégítenek bizonyos monotonitási feltételeket, illetve melyek állapotterére de niálhatunk valamilyen rendezést. Az ötlet bemutatásához tekintsük a következ® példát: Legyen
k
rögzített, az állapottér
I = 1, . . . , k ,
pij i,j∈I =
1 , 2 1 , 2
0,
32
az átmenetmátrix pedig:
ha
i, j = 1
ha
|i − j| = 1
egyébként.
Szemléletesen, ha az állapottér elemeit
1,2,. . .,k sorban vesszük, akkor a lánc minden
1 1 - valószín¶séggel, 2 2 1 illetve ha a legnagyobb vagy a legkisebb index¶ állapotban vagyunk, akkor valószí2 1 valószín¶séggel az eggyel kisebb illetve nagyobb index¶ n¶séggel helyben marad és 2 lépésben eggyel nagyobb vagy alacsonyabb index¶ állapotba lép
állapotba lép. Könnyen láthatjuk, hogy a lánc stacionárius eloszlása:
πi =
1 k
i = 1, . . . , k.
Ezek alapján legyen az átmenetszabály a következ®:
( φ(1, u) = ( φ(k, u) = valamint
i = 2,. . .,(k − 1)
1
ha
2
ha
u ∈ 0, 21 u ∈ 21 , 1 ,
k−1
ha
k
ha
u ∈ 0, 12 u ∈ 12 , 1 ,
esetén:
( φ(i, u) =
i−1
ha
i+1
ha
u ∈ 0, 21 u ∈ 12 , 1 .
Ha alkalmazzuk a Propp-Wilson algoritmust ezzel az átmenetszabállyal, lehetséges állapottal, egy jellemz® output lehet a 6.4.1.
k = 6
ábrán szerepl® illusztráció.
Ezen az ábrán szembet¶n®, hogy nincsenek egymást keresztez® átmenetek, azaz ha a
t = i-edik
id®pillanatban az egyik állapotból felfelé illetve lefelé lépünk, akkor ez
így lesz az összes többi állapot esetében is, ha azok nem a maximális vagy minimális index¶ek.
Utóbbi esetben a lánc helyben marad.
Ennek oka, hogy a fent de niált
átmenetszabály meg®rzi az állapottér rendezését: minden
{1, . . . , k}
esetén, ha
i ≤ j,
u ∈ [0, 1]
és minden
i, j ∈
akkor
φ(i, u) ≤ φ(j, u). Ez egyszer¶ következménye annak, hogy minden lánc esetében ugyanazt az használjuk.
Ha
u ∈ [0, 21 )
φ(i, u) ≤ φ(j, u).
Az
u ∈
és
i ≤ j,
akkor
i ≥ φ(i, u)
[ 21 , 1] eset hasonlóan.
valamint
ezáltal
Egyenl®ség abban az esetben állhat
fent, ha az egyik állapot maximális vagy minimális index¶.
33
j ≥ φ(j, u),
u-t
6.4.1. ábra: A Propp-Wilson algoritmus egy kimenete
k = 5 esetben, a fenti átmenetszabállyal. Az egyesüléshez a
t = −12 id®pillanatban kellett indítanunk a láncot, a lineáris rendezés szerinti legnagyobb index¶ kiindulási állapot i = 4, a legkisebb pedig j = 1, az innen induló láncokat jelöli a vastag vonal.
Ebb®l adódóan azon láncok, melyek valamely nak, mindig az
1
és
k -adik
neve is. Emiatt ha az
1
és
i ∈ {2, . . . , (k − 1)}
állapotból indul-
állapotból indulók között maradnak, innen ered az eljárás
k -adik
index¶ állapotból induló láncok egyesülnek, minden
köztes index¶nek is egyesülnie kell velük, ezáltal zéséhez elég csak kett®t nyomon követni.
k
darab lánc egyesülésének ellen®r-
A 6.4.1.
ábrán illusztrált esetben persze
nem jelentene gondot minden lánc nyomon követése, de nagy
k
esetén ez egy jelent®s
egyszer¶sítés a számítási kapacitás korlátossága miatt. Az ábrán továbbá szembet¶n®, hogy az egyesülés a legkisebb index¶ állapotban történt meg. Az átmenetszabálynak köszönhet®en az egyesülés mindig a maximális vagy minimális index¶ állapotban történik.
A függelékben található R program segítségével bemutatunk pár szimulációs
eredményt erre a feladatra:
Ha
k=5
állapotunk, akkor
106 -on
szimuláció elvégzése után a kisorsolásának szá-
mát mutatja a következ® táblázat:
k: állapot
1
2
3
4
5
darabszám
200233
200395
199836
199717
199819
A számokból jól látszik, hogy az eloszlás közelít®leg egyenletes.
34
A következ® táblázatban különböz®
k állapotszámokra 50-50 szimuláció eredménye-
inek átlagait láthatjuk: k
3
4
5
6
10
20
Bolyongás hosszának átlaga
2.94
5.34
10.5
14.06
40.2
216.98
Egyesülés idejének átlaga
2
3.64
6.04
8.64
27.36
124.2
Egyesülés utáni bolyongási hossz átlaga
0.94
1.7
4.46
5.42
12.84
92.78
A táblázat eredményei alapján a bolyongás hosszának átlaga nem lineáris függvénye az állapotszámnak. Nyilvánvaló, hogy
k
állapot esetén legalább
(k − 1)
lépést kell
tennünk, hogy a fels® és az alsó láncok találkozzanak. A két szám kapcsolatából viszont arra következtethetünk, hogy az algoritmus futási idejének csökkentése céljából érdemes nagy
k
T
esetén nagyobb kezd®
értéket választanunk, így elhagyva azokat a
lépésszámokat, melyek esetén elhanyagolható a találkozás valószín¶sége, csökkentjük a végigszámolt esetek számát. A 6.1.-es ábra az egyesülési id®ket mutatja a bolyongások teljes hosszának függvényében. Itt a
k = 4, 7, 10, 20
állapotszámú eseteket láthatjuk 50-50 szimuláció után.
Jól látszik, hogy a bolyongások hossza gyorsabban növekszik, mint az egyesülési id®, ha
k -t
növeljük.
Az eddigiek alapján a módszer általános leírásához legyen adott egy
≤
részben rendezettség, valamint a
tonitási feltétel, azaz
x≤y
fel továbbá, hogy az
I
esetén
állapottér
φ
I
az állapottér, melyen
átmenetszabályra teljesüljön a mono-
φ(x, U0 ) ≤ φ(y, U0 ), minden x, y ∈ I
x
elemeire:
0 ≤ x ≤ 1.
esetén. Tegyük
Legyen
Φtt21 (x, u) = φt2 −1 (φt2 −2 (. . . (φt1 (x, ut1 ), ut1 +1 ), . . .), ut2 −2 ), ut2 −1 ) ahol
u = (. . .,u−1 , u0 ).
Ha
u−T ,u−T +1 ,. . .,u−2 ,u−1 -re
teljesül, hogy
Φ0−T (0, u) =
Φ0−T (1, u), akkor a monotonitási feltétel biztosítja, hogy Φ0−T (x, u) = Φ0−T (y, u), minden x, y ∈ I T -vel, után a
T,
esetén. Ezáltal a legkisebb
melyre
melyre
Φ0−T (0, u) = Φ0−T (1, u).
T = 1,2,4,8,. . .
Φ0−T (·, u)
Jelöljük ezt a
konstans, egyenl® a legkisebb
T -t T ∗ -gal.
értékeket, amíg nem találunk egy
k -t,
Próbáljuk egymás
melyre
T = 2k
esetén
Φ0−T (0, u) = Φ0−T (1, u). Ezáltal a szimuláció lépésszáma 2(1 + 2 + 4 + . . . + 2k ) < 2k+2 , ahol az egyenl®tlenség elején szerepl® 2-es szorzó a két lánc párhuzamos futtatása miatt jelenik meg (az egyiket a optimális, ugyanis
T
∗
0-ból,
-nak legalább
k−1
2
a másikat az
1-b®l
indítjuk).
Ez közelít®leg
-nek kell lennie, különben nem kellett volna
35
T =
(a) k=4 állapot
(c)
k =10
(b)
állapot
(d)
k =7
k =20
állapot
6.1. ábra. Az egyesülés ideje a bolyongások hosszának függvényében
36
2k
lépést tennünk, ezáltal a szükséges szimulációs lépések száma, hogy megmutassuk
Φ0−T ∗ (0, u) = Φ0−T ∗ (1, u)
nagyobb mint
2 · 2k−1 = 2k .
Az eljárás összegzése algoritmikus formában:
T
←
1
REPEAT
upper ← 1 lower ← 0 FOR
t = −T
to
−1
upper ← φt (upper, ut ) lower ← φt (lower, ut ) T ← 2T UNTIL
upper = lower
RETURN
upper
6.2. Megjegyzés. Mint a korábbi, általános CFTP esetében, itt is fontos, hogy ugyan-
azokat az u-kat használjuk, és ne generáljunk a REPEAT során újakat.
Mint láthattuk, a sandwiching módszer a gyakorlatban jól alkalmazhatóvá teszi a CFTP algoritmust, de sajnos olyan Markov láncok eseten nem alkalmazható, melyek átmenetszabályára nem teljeseül a monotonitási feltétel, vagy nem adható meg egy részben rendezettség maximális és minimális állapotokkal az állapottéren.
37
7. fejezet Összefoglalás
A dolgozat elkészítésével az volt a célom, hogy megismerjem az addig számomra kevéssé ismert Markov láncokat és az azokat felhasználó MCMC módszereket. Az els® fejezetben bemutatásra kerültek a Markov láncok, mely jórészt az
[1]-en
alapult. Ezután az
algoritmus család id®beli fejl®dését követve bemutattam a talán leginkább elterjedet eljárásokat, kezdve a rengeteg helyen alkalmazott Metropolis-Hastings algoritmussal és a Gibbs mintavételez®vel. Ezen fejezetek alapját az
[5],[6],[7] szolgáltak.
Némi kitekin-
tésekt®l eltekintve, a dolgozat kereteinek fényében szinte mindenhol diszkrét eseteket tárgyaltam.
Ezt követte a Slice sampling általános eljárásának bemutatása illetve a
megvalósítás különböz® lehet®ségeinek tárgyalása. A fejezet alapját a
[16]
képezi. Vé-
gül pedig bemutatásra került a Perfect sampling ötlete, egy általánosan használható eljárása és egy, a gyakorlati alkalmazást is lehet® tév® kiegészítése, a sandwiching módszer. Ezt a módszert egy megértést el®segít® példán keresztül mutattam be, illetve a Perfect Samplingre jellemz® tulajdonságok vizsgálatára egy R programot is készítettem.
Ez a program kiterjeszthet® további eloszlások mintavételezésére is.
alapját a
[9],[10],[12]
és
[14]
A fejezet
adta. A Perfect samplinggel kapcsolatban még manapság
is számos nyitott probléma van ezért sok kutatás tárgyát képezi.
38
8. fejezet Függelék
8.1.
R program a sandwiching módszer bemutatásához
Az alábbi programkód a dolgozatban szerepl® eredmények kiszámítását végzik el. Ha azok csak közül 1-1-re vagyunk kíváncsiak, a megfelel® részeket kommentként átírhatjuk. A sandw függvény for ciklusba ágyazása miatt, ha csupán az egészet bemásoljuk az R-be, nem jelenik meg a DF és a DF2 adatsor, mely az állapotok változásának folyamatát írja ki a konzolra.
setwd('C:/R') ### Creating a uniform distribution on {0,1} distr <- function() { randomNumber <- runif(1, min=0, max=1) if (randomNumber >= 0.5) randomNumber = 1 else if (randomNumber < 0.5) randomNumber = -1 } nb <- 50 # set the number of simulations 39
c <- c(1:nb) d <- c(1:nb) e <- c(1:nb) for (i in 1:nb) { upd <- distr() # set upd's beginning value k <- 20 # set k as the number of possible states ### Building Sandwiching function sandw <- function (upd) { x <- seq(1:k) # index vector mtx <- matrix(0,k,1) # empty state matrix DF = data.frame(x) mav <- k # state with maximal index miv <- 1 # state with minimal index y <- seq(1:k) mtx2 <- matrix(0,k,1) DF2 = data.frame(y) ma <- k mi <- 1 for (i in 1:length(upd)) { mtx[mav,1] <- 0 # set previous element with index mav to 0 mtx[miv,1] <- 0 # set previous element with index miv to 0 mav <- min(max(mav+upd[i],1),k) # update mav miv <- min(max(miv+upd[i],1),k) # update miv
40
mtx[mav,1] <- 1 # update the vector element mtx[miv,1] <- 1 # update the vector element DF <- data.frame(mtx,DF) if (ma!=mi) # do until reach the state of coalescence { mtx2[ma,1] <- 0 mtx2[mi,1] <- 0 ma <- min(max(ma+upd[i],1),k) mi <- min(max(mi+upd[i],1),k) mtx2[ma,1] <- 1 mtx2[mi,1] <- 1 DF2 <- data.frame(mtx2,DF2) } } if (mav!=miv) return(NA) else list(mav,DF,DF2,ncol(DF)-1,ncol(DF2)-1) # print random walks and numbers of them } while (is.na(sandw(upd))) upd <- c(distr(), upd) # do until reaching the state of coalescence sandw(upd)
# listing
# create a vector with maximum values of the simulations c[i] <- sandw(upd)[1] # create a vector with total lengths of random walks d[i] <- sandw(upd)[4] # create a vector with number of steps by reaching coalescence e[i] <- sandw(upd)[5] 41
} d <- unlist(d) e <- unlist(e) table(unlist(c)) # summary of draws tb <- data.frame(d,e,d-e) mean1 <- sum(tb[1])/nb # mean: total lengths of random walks mean2 <- sum(tb[2])/nb # mean: number of steps by reaching coalescence mean3 <- sum(tb[3])/nb # mean: number of steps after reaching coalescence mean1; mean2; mean3 jpeg(filename="plot.jpg") matplot(d,e,type="p",lwd=2,col=3,xlab="lengths of random walks",ylab= "number of steps by reaching coalescence") dev.off()
42
Irodalomjegyzék
[1] http://www.stat.u .edu/ casella/Papers/MCMCHistory.pdf [2] http://www.cs.elte.hu/ villo/ml/ML.pdf [3] I. M. Szobol: A Monte-Carlo mószerek alapjai
M¶szaki Könyvkiadó, Budapest,
1981 [4] Gilks W.R., Richardson S., Spiegelhalter D.J.: Carlo in Practice
Chapman & Hall/CRC, 1996.
Markov Chain Monte
[5] http://web.mit.edu/ wingated/www/introductions/mcmc-gibbs-intro.pdf [6] Siddhartha Chib,
Edward
Hastings Algorithm. In:
Greenberg:
Understanding the Metropolis-
The American Statistician, Vol. 49, No. 4. (Nov., 1995),
pp. 327-335 [7] George Casella, Edward I. George:
Explaining the Gibbs Sampler. In:
The American Statistician, Vol. 46, No. 3 (Aug., 1992), pp. 167-174.
[8] http://www.stat.lsa.umich.edu/ kshedden/Courses/Stat606/Notes/mcmc.pdf [9] George Casella, Michael Lavine, Christian P. Robert: Explaining the Perfect Sampler. In:
The American Statistician, Vol. 55, No. 4. (Nov., 2001), pp.
299-305. [10] Olle Haggstrom: Finite Markov Chains and Algorithmic Applications
ridge University Press, 2002.
[11] David B. Wilson: Exact Sampling with Markov Chains
of Technology, 1996.
43
Camb-
Massachusetts Institute
[12] Wilfrid S. Kendall: Notes on Perfect Simulation
University of Warwick, 2004.
[13] Solymosi Norbert: Bevezetés az R-nyelv és környezet használatába, 2005. [14] James J. Propp, David B. Wilson: Exact Sampling with Coupled Markov Chains and Applications to Statistical Mechanics.
nology, 1996.
Massachusetts Institute of Tech-
[15] László Lovász, Peter Winkler: Exact Mixing in an Unknown Markov Chain.
The Electronic Journal of Combinatorics, Vol. 2, R15 (1995)
[16] Radford M. Neal: Slice sampling. In: (Jun., 2003), pp. 705-741.
44
The Annals of Statistics, Vol. 31, No. 3
9. fejezet Köszönetnyilvánítás
Ezúton szeretném megköszönni témavezet®mnek, Pröhle Tamásnak a dolgozat elkészítésében nyújtott segítségét, folyamatos útmutatását és hogy bevezetett az R program használatába. Köszönöm évfolyamtársamnak, Csillagvári Dánielnek az ábrák elkészítésében nyújott segítségét.
45
Nyilatkozat
A szakdolgozat szerz®jeként fegyelmi felel®sségem tudatában kijelentem, hogy a dolgozatom önálló munkám eredménye, saját szellemi termékem, abban a hivatkozások és idézések standard szabályait következetesen alkalmaztam, mások által írt részeket a megfelel® idézés nélkül nem használtam fel.
46
