Monte Carlo Markov Chains voor Bayesiaanse statistiek

Bachelorscriptie

Monte Carlo Markov Chains voor Bayesiaanse statistiek Jochem Braakman, 19 juli 2012 Begeleiding: dr. B.J.K. Kleijn

1 0.9 0.8 0.7

µ

0.6 0.5 0.4 0.3 0.2 0.1 0

0

1

2

3

4

5 λ

6

7

8

9

10

KdV Instituut voor wiskunde Faculteit der Natuurwetenschappen, Wiskunde en Informatica Universiteit van Amsterdam

Samenvatting In deze scriptie nemen wij een kijkje in de wereld van Monte Carlo Markov Chains. Monte Carlo Markov Chains zijn een klasse van algoritmes die ons in staat stellen kansverdelingen te simuleren. De simulatie wordt uitgevoerd met behulp van een computer. In de Bayesiaanse statistiek wil men de zogenaamde posterior verdeling bepalen. Omdat de posterior in de meeste situaties niet in een gesloten vorm beschikbaar is, wil men deze discreet benaderen. Dit kan met behulp van Monte Carlo Markov Chains. De scriptie is als volgt opgebouwd. In het eerste hoofdstuk gaan we naar de basis van simulatietechnieken en gaan we in op Markov Chains. In het tweede hoofdstuk gaan we in op de beginselen van de Bayesiaanse statistiek om de posterior verdeling af te leiden. Daarna leggen wij uit waarom de posterior zo belangrijk is voor een Bayesiaanse statisticus. In het derde hoofdstuk combineren wij de theorie van Hoofdstuk 1 en Hoofdstuk 2 om vervolgens de Monte Carlo Markov Chains te introduceren. We bespreken daarbij twee methodes: Het symmetrische Metropolis Hastings algoritme en de Gibbs Sampler. Aan het eind van het hoofdstuk voeren we simulaties van bekende en onbekende kansverdelingen met behulp van de reeds genoemde algoritmes uit.

Gegevens Titel: Monte Carlo Markov Chains voor Bayesiaanse statistiek Auteur: Jochem Braakman, [email protected], 6035868 Begeleider: dr. B.J.K. Kleijn Tweede beoordelaar: dr. Bert van Es Einddatum: 19 juli 2012 Korteweg de Vries Instituut voor Wiskunde Universiteit van Amsterdam Science Park 904, 1098 XH Amsterdam http://www.science.uva.nl/math

Inhoudsopgave 1 Inleiding 1.1 Monte Carlo integratie . . . . . 1.2 Inversie en Verwerpingsmethode 1.2.1 Verwerpingsmethode . . 1.3 Markov Chains . . . . . . . . . 1.4 Convergentie Markov Chains . .

. . . . .

. . . . .

. . . . .

. . . . .

. . . . .

. . . . .

. . . . .

. . . . .

. . . . .

. . . . .

. . . . .

. . . . .

. . . . .

. . . . .

. . . . .

. . . . .

. . . . .

2 2 4 6 9 13

2 Bayesiaanse statistiek 15 2.1 Prior en posterior verdelingen . . . . . . . . . . . . . . . . . . 15 2.2 Inferentiële aspecten van de posterior . . . . . . . . . . . . . . 19 3 Monte Carlo Markov Chains 23 3.1 Symmetrisch Metropolis Hastings algoritme . . . . . . . . . . 23 3.2 Gibbs Sampler . . . . . . . . . . . . . . . . . . . . . . . . . . 25 3.3 Voorbeelden . . . . . . . . . . . . . . . . . . . . . . . . . . . . 28 4 Conclusie en discussie 38 4.1 Samenvattend . . . . . . . . . . . . . . . . . . . . . . . . . . . 38 4.2 Resterende vragen . . . . . . . . . . . . . . . . . . . . . . . . 40 Populaire samenvatting

41

Bibliografie

43

1

Hoofdstuk 1 Inleiding In deze scriptie zullen we ons richten op Monte Carlo Markov Chains en hun toepassingen in de Bayesiaanse statistiek. De Bayesiaanse statisticus is ge¨ınteresseerd in de posterior verdeling. De posterior verdeling stelt de Bayesiaan in staat inferentiële aspecten van de statistiek zoals bijvoorbeeld schatten, toetsen, het vinden van betrouwbaarheidsverzamelingen te kunnen doen. In het algemeen is het niet zo eenvoudig de posterior in een gesloten vorm te vinden en daarom is ons doel van de scriptie de posterior numeriek te benaderen. Hier bieden Monte Carlo Markov Chains uitkomst. Voor verdere informatie over de posterior en het belang daarvan verwijzen wij naar hoofdstuk 2. Voor de Monte Carlo Markov Chains verwijzen wij naar hoofdstuk 3. In dit inleidende hoofdstuk houden we ons bezig met het simuleren van kansverdelingen en de theorie van Markov Chains.

1.1

Monte Carlo integratie

R Zij een ruimte X gegeven. Stel wij willen een integraal X g(x) dx berekenen. Als het mogelijk is g te factoriseren in een functie f (x) en een kansdichtheid π(x) dan geldt: Z Z g(x) dx = f (x)π(x) dx = Eπ (f (x)). X

X

De wet van de grote aantallen implicieert, dat als men een groot aantal X1 , ..., Xn uit de kansdichtheid π kan simuleren, dat de verwachting benaderd kan worden door het steekproefgemiddelde. n

1X f (Xi ) −→ Eπ (f (X)), n i=1 2

(π −bijna zeker ).

De Monte Carlo integreermethode stelt ons in staat integralen van bovenstaande vorm te benaderen. De vraag die direct volgt aangaande de nauwkeurigheid van deze benadering, wordt beantwoord door de centrale limietstelling. Als men i.i.d. X1 , ..., Xn stochasten trekt met marginale kansdichtheid π en als deze eindige variantie σ 2 hebben, geldt dat, √ n

n

1X Xi − EX1 n i=1

! N (0, σ 2 ).

We concluderen dat het√verschil tussen de approximatie en de verwachting van X1 met grofweg 1/ n naar 0 gaat. Voorbeeld 1.1. We willen de variantie met schatten voor de stochast X die de Student t-verdeling heeft met een onbekend aantal ν > 2 vrijheidsgraden. Deze heeft dichtheid ) 1 Γ( ν+1 2 √ f (x|ν) = ν Γ( 2 ) νπ

x2 1+ ν

−( ν+1 2 ) .

De verwachting EX is bekend en wordt gegeven door 0, zodat de variantie van X gelijk is aan EX 2 . Met de Monte Carlo benadering, Z

n

∞

1X 2 x f (x|ν) dx = Ef (x|ν) (X ) ≈ X , EX = n i=1 i −∞ 2

2

2

waarbij wij X1 , ..., Xn willen trekken uit f (.|ν). De resterende vraag is hoe wij punten kunnen simuleren zo dat deze bij benadering volgens de Student t-verdeling verdeeld zijn. Dat antwoord zullen wij in de populaire samenvatting geven. Daar zullen we ook met behulp van zo genaamde Monte Carlo Markov Chains de varianties van de Student t-verdeling voor verschillende vrijheidsgraden ν schatten. Opmerking 1.2. Voor dat wij ons met algoritmes bezig houden, beginnen wij bij de basis van simulatie. Alles wat wij hier bij nodig hebben is een random number generator. Een random number generator genereert getallen tussen 0 en 1 die uniform U [0, 1] verdeeld zijn. Omdat computerpaketten zoals MATLAB over random number generators beschikken houden we ons in deze scriptie niet met de vraag bezig, hoe een dergelijke random number generator functioneert en beschouwen dit als een gegeven.

3

1.2

Inversie en Verwerpingsmethode

We zullen nu twee simulatiemethodes geven om uit een bepaalde distributie willekeurige getallen te trekken. Gegeven een random number generator die uniform getallen u tussen 0 en 1 genereert, geven wij de inversiemethode voor continue stochasten. Zij F de verdelingsfunctie waaruit we een trekking willen simuleren. Als F een strict stijgende, continue verdelingsfunctie heeft, definiëer dan de stochast X door: X ∼ F −1 (U ). Stelling 1.3. Gegeven een verdelingsfunctie F . Zij (U1 , ..., Un ) ∼ U [0, 1]n en definiëer Xi = F −1 (Ui ) voor 1 ≤ i ≤ n. Dan is (X1 , ..., Xn ) ∼ F n . Bewijs. Aangezien voor alle A = A1 × ... × An geldt dat, P(X1 ∈ A1 , ..., Xn ∈ An ) = P(U1 ∈ F (A1 ), ..., Un ∈ F (An )) n n Y Y = P(Ui ∈ F (Ai )) = P(Xi ∈ Ai ), i=1

i=1

vormen de X1 , ..., Xn een onafhankelijke steekproef. Uit het feit dat P(Xi ≤ x) = P(F −1 (Ui ) ≤ x) = P(Ui ≤ F (x)) = F (x), volgt daarnaast, dat Xi marginaal verdeeld is volgens F voor iedere 1 ≤ i ≤ n. Een i.i.d. steekproef U1 , ..., Un ∼ U [0, 1]n kan met behulp van deze stelling worden getransformeerd naar een i.i.d. steekproef X1 , ..., Xn ∼ F n . Voorbeeld 1.4. Beschouw voor zekere λ > 0 de exponentiële verdeling met dichtheid fλ (x) = λe−λx 1{x≥0} . De bijbehorende verdelingsfunctie is Fλ (x) = (1 − eλx )1{x≥0} met inverse Fλ−1 (u) = − log(1−u) voor u ∈ (0, 1). Toepassing λ van stelling 1.3 ter verkrijging van een gesimuleerde steekproef X1 , ..., Xn leidt tot het algoritme van figuur 1.1. Daadwerkelijke simulatie van een i.i.d. steekproef ter groote N = 20000 leidt tot het histogram in figuur 1.2.

4

1. Trek een getal u uit U [0, 1]

f u n c t i o n [X] = i n v e r s i e e x p o n e n t i e e l ( lambda ,N) f o r i =1:N

2. Bereken F −1 (u) = − log(1−u) λ 3. Herhaal tot een bovengrens N.

U=rand ; X( i )=−l o g (1−U) / lambda ; end

end

Figuur 1.1: Links zien wij een beschrijving van het inversiealgoritme. Rechts implementeren wij deze beschrijving in een MATLAB code.

4000 3500 3000 2500 2000 1500 1000 500 0

0

2

4

6

8

10

12

Figuur 1.2: Histogram met N = 20000 trekkingen door middel van de inversiemethode voor de exponentiële verdeling met λ = 1. Opmerking 1.5. De inversiemethode werkt voor alle verdelingen waarvoor de verdelingsfunctie inverteerbaar is, indien de inverse F −1 in gesloten vorm beschikbaar is. Dit is niet het geval voor bijvoorbeeld de normale verdeling. De verdelingsfunctie van de normale verdeling is een integraal die geen 5

primitieve heeft. Dit heeft tot gevolg dat we geen expliciete vorm van de inverse kunnen vinden voor de verdelingsfunctie van de normale verdeling. Een sterker voorbeeld wordt gevormd door Bayesiaanse posteriors. In Bayesiaanse context is het probleem ook aanwezig, aangezien de posterior weliswaar een formele definitie heeft, maar helemaal niet eenvoudig toegankelijk is in praktisch (of computationele) zin. Het is dus zaak naar andere methodes te kijken. Een bekend algoritme voor trekkingen uit een verdeling is de verwerpingsmethode. Net als bij de inversie werken we met continue stochasten.

1.2.1

Verwerpingsmethode

Stel wij willen random getallen uit een verdeling F simuleren en we kunnen random getallen genereren uit een verdeling G. We noteren de dichtheden behorende bij de verdelingen f respectievelijk g. Zij M > 0 een constante zo (x) ≤ M voor alle x en zo dat g(x) 6= 0 voor alle x. We hebben dan het dat fg(x) volgende algoritme om een stochast X met dichtheid f te simuleren. 1. Genereer een waarde Y uit G en trek onafhankelijk daarvan een U ∼ U [0, 1]. 2. Als U ≤

f (Y ) M g(Y )

dan is X = Y . Wij accepteren Y .

3. Als U >

f (Y ) M g(Y )

verwerp Y, wijs niets toe aan X en begin opnieuw.

Deze procedure kunnen wij tot een vaste bovengrens N herhalen om een i.i.d. steekproef te genereren. Stelling 1.6. De stochast X die gegenereerd wordt door de verwerpingsmethode heeft dichtheid f . Bewijs. P(X ≤ x) = P(Y ≤ x|Y geaccepteerd ) (Y ) P(X ≤ x, U ≤ Mf g(Y ) f (Y ) ) = P X ≤ x U ≤ = . M g(Y ) P(U ≤ f (Y ) )) M g(Y )

6

Zij m = P(U ≤

f (Y ) )). M g(Y )

Dan is

P(X ≤ x) = =

P(X ≤ x, U ≤

P(X ≤ x, U ≤

m f (Y ) |Y M g(Y )

= y)P(Y = y)

m Rx

=

f (Y ) ) M g(Y )

f (y) g(y) dy −∞ M g(y)

m

Rx =

f (y) −∞ M

m

dy

.

Neem de limiet x → ∞, dan R∞ P(X ≤ x) = 1 =

f (y) dy , Mm

−∞

hetgeen impliceert dat m = 1/M . Opmerking 1.7. Het beste zou zijn als de acceptatiekans m = P(U ≤ f (Y ) ) = 1 is. Dit zou betekenen dat f = g bijna overal. Maar dat is een M g(Y ) tegenspraak met de aanname dat we uit de verdeling G kunnen simuleren en niet uit F . Dus M > 1. Merk op dat 0 ≤ f /(gM ) ≤ 1. Wij noemen dan f /(gM ) ook de acceptatiefunctie voor het verwerpingsalgoritme. Het algoritme is het meest effectief als M klein is, want dan is de kans op acceptatie hoog en worden tijdrovende herhalingen met hoge kans vermeden. De verwerpingsmethode stelt ons in staat om uit de normale verdeling trekkingen te doen, hetgeen van praktisch nut zal zijn in het vervolg van deze scriptie. D

Voorbeeld 1.8. Als X ∼ N (µ, σ 2 ) dan is X = σZ + µ, waarbij Z ∼ N (0, 1) de standaard normale verdeling is. Het is voldoende om uit de standaard normale verdeling te trekken. Omdat de standaard normale verdeling een symmetrische dichtheid rond 0 heeft, is het voldoende om uit |Z| trekkingen te doen en het teken van Z te trekken uit de Bernoulli( 21 ) verdeling. De x2

stochast |Z| heeft dan dichtheid f (x) = √22π e− 2 1{x≥0} . We kiezen nu een verdelingsfunctie G met dichtheid g met de eigenschap dat er een M > 1 bestaat zodat f ≤ M g en dat we makkelijk uit G kunnen simuleren. Hiertoe bedienen wij ons uit het vorige voorbeeld, namelijk met G ∼ exp(1). Volgens (x) het verwerpingsalgoritme moeten we h(x) = Mf g(x) kiezen, waarbij we M p f (x) x−x2 /2 moeten berekenen. Nu is g(x) = e 2/π. Het maximum M van de p functie ligt bij x = 1. Dan is M = f (1)/g(1) = 2e/π. We krijgen dan 2 h(x) = e−(x−1) /2 . Om uiteindelijk uit Z te trekken definiëren we een stochast 7

U ∼ U [0, 1] en accepteren wij Z = |Z| als U ≤ 0.5 en zetten we Z = −|Z| als U > 0.5. Dit doen we omdat de kansmassa van Z symmetrisch rond 0 ligt. Het algoritme om uit Z te trekken ziet er dan als volgt uit: 1. Trek Y uit exp(1). 2. Genereer U uit U [0, 1]. 3. Accepteer |Z| = Y als U ≤ e−(Y −1)

2 /2

ga terug naar de eerste stap.

4. Genereer een nieuwe U uit U [0, 1], zet Z = |Z| als U ≤ 0.5 en zet Z = −|Z|. Opmerking 1.9. We kunnen dit algoritme nog versimpelen. Omdat U ≤ 2 e−(Y −1) /2 is dit equivalent met − log(U ) ≥ (Y − 1)2 /2. Nu is volgens de inversiemethode − log(U ) ∼ exp(1). Het resultaat is een eenvoudiger algoritme:

1. Trek Y en U onafhankelijk van elkaar uit exp(1)

f u n c t i o n [X] = r e j e c t i o n m e t h o d n o r m a l (N, sigma , mu) i =1; w h i l e i <=N Y=i n v e r s i e e x p o n e n t i e e l ( 1 , 1 ) ; U=i n v e r s i e e x p o n e n t i e e l ( 1 , 1 ) ; i f U>=0.5∗(Y−1) ˆ2 Z ( i )=Y; i=i +1; end end

2. accepteer |Z| = Y als U ≥ (Y − 1)2 /2 ga anders terug naar stap 1 3. Genereer een nieuwe U uit U [0, 1], zet Z = |Z| als U ≤ 0.5 en zet Z = −|Z| elders.

f o r i =1: l e n g t h ( Z ) U=rand ; i f U<=0.5 Z ( i )=Z ( i ) ; else Z ( i )=−Z ( i ) ; end end X=sigma ∗Z+mu ; h i s t (X, 3 0 ) end

Figuur 1.3: Links zien we een beschrijving van het algoritme om met behulp van de verwerpingsmethode uit een normale verdeling N (µ, σ 2 ) te trekken en rechts zien wij de implementatie van de beschrijving in een MATLAB code. 8

1800 1600 1400 1200 1000 800 600 400 200 0 −5

−4

−3

−2

−1

0

1

2

3

4

5

Figuur 1.4: Histogram voor N = 10000 trekkingen uit de standard normale verdeling N (0, 1) met behulp van de verwerpingsmethode. We hebben gezien dat de basis voor simulatie een random number generator is. Vanuit daar kan men de exponentiële verdeling simuleren. Deze verdeling hebben wij op zijn beurt gebruikt om met de verwerpingsmethode de normale verdeling te simuleren. De vraag of men daarmee uit de posterior kan simuleren moet met nee beantwoord worden. Toch vormen deze twee onderdelen een belangrijke basis voor het nog komende. We zullen onze theorie over simulatie moeten uitbreiden om in staat te zijn uit de posteriors trekkingen te doen.

1.3

Markov Chains

In hoofdstuk 3 passen we Monte Carlo Markov Chains toe op posteriors. Zoals de naam Monte Carlo Markov Chains al suggereert moeten we iets te weten komen over Markov Chains. Dit doen we door middel van definities en stellingen om het uiteindelijke doel, de convergentie naar een evenwichtsdistributie, te bereiken. (Stelling 1.27). Omwille van de lengte en de leesbaarheid van deze sectie zullen we de stellingen niet in detail bewijzen en meestal naar het boek Markov Chains van James R. Norris verwijzen [5]. De 9

resultaten die we krijgen stellen ons in staat door de combinatie van Markov Chains en de simulatiemethoden uit hoofdstuk 1 een Monte Carlo Markov Chain te construeren, met die we uit posterior verdelingen kunnen trekken. Zie verder hoofdstuk 2 en 3. Definitie 1.10. Zij I een aftelbare verzameling. We noemen i ∈ I een toestand en I de toestandsruimte. We noemen P λ = (λi : i ∈ I) een maat als geldt dat 0 ≤ λi < ∞ voor alle i. Indien i∈I λi = 1, dan wordt λ een distributie genoemd. Definitie 1.11. We noemen een matrix P = (pij : i, jP ∈ I) de transitiematrix van de Markov Chain als pij ≥ 0 en voor alle j ∈ I, i∈I pij = 1. Definitie 1.12. We zeggen dat (Xn )n≥0 een Markov Chain is met begindistributie λ en transitiematrix P , als er aan de volgende twee voorwaarden voldaan wordt: (i) X0 is verdeeld volgens distributie λ. (ii) voor n ≥ 0, gegeven Xn = i, heeft Xn+1 de distributie (pij : j ∈ I), i.e. Xn+1 is onafhankelijk van X0 , ..., Xn . Met gebruik van de notatie P voor de kansmaat geassociëerd met (Xn ) kunnen we deze definiërende eigenschappen voor toestanden i1 , ..., in+1 ∈ I herformuleren als volgt: (i) P(X0 = i0 ) = λi0 , (ii) P(Xn+1 = in+1 |X0 = i0 , ...., Xn = in ) = P(Xn+1 = in+1 |Xn = in ) = pin in+1 . We zeggen dan ook wel dat het proces (Xn )n≥0 een Markov-proces is beginnend bij λ, met transitiematrix P . We noteren dit ook wel met Markov(λ, P ). Met de gebruikelijke definitie van de Kronecker delta δij = 1{i = j} P beschouwen wij nu (λP )j = i∈I λi pij . Stelling 1.13. (Markov eigenschap) Laat (Xn )n≥0 Markov(λ, P ) zijn. Gegeven Xm = i, is (Xm+n )n≥0 Markov(δi , P ) en onafhankelijk van X0 , ..., Xm . Bewijs. Het bewijs is te vinden in het boek ”Markov Chains” (Norris (1998) [5]). (n)

We voeren de volgende notatie in: pij is de (i, j)-de entry in de matrix (n) P n , de n-staps transitiematrix. We noemen pij de n-staps overgangskans tussen i en j. Verder als λi > 0 noteren we P(A|X0 = i) met Pi (A). Door 10

de Markov eigenschap op tijdstip m = 0, onder Pi is (Xn )n≥0 Markov(δi , P ). Dus (Xn )n≥0 onder Pi hangt niet af van λ. Stelling 1.14. Zij (Xn )n≥0 Markov(λ, P ). Dan geldt voor alle m, n ≥ 0 dat i) P(Xn = j) = (λP n )j , (n)

ii) Pi (Xn = j) = P(Xn+m = j|Xm = i) = pij . Bewijs. Er geldt voor i) door over alle toestanden i1 , ..., in−1 te sommeren dat X X P(Xn = j) = ... P(X0 = i1 , ..., Xn−1 = in−1 , Xn = j) i1 ∈I

=

X i1 ∈I

in−1 ∈I

...

X

λi1 pi1 i2 ...pin−1j = (λP n )j .

in−1 ∈I

Voor ii) gebruiken wij de Markov eigenschap uit de 1.13. Die zegt gegeven Xm = i, dat (Xm+n )n≥0 Markov(δi , P ) is. We kiezen om ii) te krijgen λ = δi in i). (n)

Wij zijn in het vervolg ge¨ınteresseerd in limn→∞ pij , de overgangskans na ∞ aantal stappen. Definitie 1.15. We zeggen dat voor i, j ∈ I dat i naar j leidt, als Pi (∃n ≥ 0 | Xn = j) > 0. We schrijven dan ook i → j. Als er nu ook geldt dat j → i, dan zeggen we dat i en j met elkaar communiceren. Wij noteren dit met i ↔ j. Stelling 1.16. Voor de toestanden i, j ∈ I zijn volgende uitspraken equivalent: (i) i → j. (ii) Er bestaat een n en i1 , ..., in zdd: pi1 i2 pi2 i3 ...pin−1 in > 0 voor toestanden i1 , ..., in ∈ I met i1 = i, in = j. (n)

(iii) pij > 0 voor een zekere n ≥ 0. Bewijs. Voor het bewijs zie ook ”Markov Chains” (Norris (1998) [5]).

11

Wij kunnen aan de hand van communicerende toestanden een equivalentierelatie ∼I geven op de toestandsruimte I door te stellen dat i ∼I j d.e.s.d.a. i ↔ j. Door n = 0 te kiezen is meteen duidelijk dat voor een i ∈ I geldt dat i ∼I i. Als geldt dat i ∼I j dan geldt dat i ↔ j hetgeen equivalent is met j ↔ i hetgeen equivalent is met j ∼I i. Dit bewijst symmetrie. Als i ∼I j en j ∼I k, betekent dit dat i ↔ j en j ↔ k. Dit betekent in het bijzonder dat i → j en j → k. Met behulp van stelling 1.16 is dit equivalent met pij > 0 en pjk > 0, maar dan is pik = pij pjk > 0 hetgeen betekent dat i → k. Op een zelfde manier kunnen wij aantonen dat k → i, hetgeen de transitiviteit bewijst. Dus ∼I is een equivalentierelatie op I die partitioneert in communicerende klassen. Definitie 1.17. Een Markov Chain (Xn )n≥0 met transitiematrix P heet irreducibel als de toestandsruimte I onder de equivalentierelate ∼I één enkele communicerende klasse is. Definitie 1.18. Een stochast T : Ω → {0, 1, 2, ...} ∪ {∞} is een stoptijd als voor iedere n ≥ 0, het evenement {T = n} alléén afhangt van X0 , ..., Xn voor n = 0, 1, 2, ... . Stelling 1.19. (De sterke Markov eigenschap.) Zij (Xn )n≥0 een Markov Chain met transitiematrix P en zij T een stoptijd. Gegeven T < ∞ en XT = i, is (XT +n )n≥0 een Markov Chain met begindistributie δi en overgangsmatrix P en is deze onafhankelijk van X0 , ..., XT . Bewijs. Het bewijs is te vinden in Norris (1998) [5]. Definitie 1.20. Een toestand i ∈ I van een Markov Chain (Xn )n≥0 is recurrent als Pi (Xn = i voor oneindig veel n) = 1. Een toestand i ∈ I heet transiënt als Pi (Xn = i voor oneindig veel n) = 0. Definitie 1.21. De first passage time naar een toestand i is de stochast Ti die gedefiniëerd wordt door Ti (ω) = inf{n ≥ 1 : Xn (ω) = i}, met inf ∅ = ∞ Stelling 1.22. Stel P is irreducibel en recurrent, dan is P(Tj < ∞) = 1. Bewijs. Voor het bewijs zie Norris (1998) [5]. 12

Ook hier is de stelling alleen van belang voor het nog komende resultaat over convergentie van Markov Chains 1.27. Definitie 1.23. Als een toestand i recurrent is en daarbij de verwachte terugkeer tijd mi = Ei (Ti ) eindig is, dan noemen wij i ook wel positief recurrent. Definitie 1.24. We noemen een maat π invariant, als geldt dat πP = π. Als bovendien geldt dat

P

i∈I

πi = 1, noemen wij π een invariante distributie.

Nu volgt weer een stelling die nodig is voor een belangrijk eindresultaat voor de convergentie van Markov Chains. Stelling 1.25. Zij P irreducibel. Dan zijn de volgende uitspraken equivalent: (i) iedere toestand is positief recurrent; (ii) een toestand i is positief recurrent; (iii) P heeft een invariante maat die we π noemen. In dat geval is mi = 1/πi voor alle i. Bewijs. Ook hier verwijzen we naar Norris (1998) [5].

1.4

Convergentie Markov Chains

Definitie 1.26. Een toestand i ∈ I heet aperiodiek als er een N ∈ N bestaat (n) zodanig dat voor alle n ≥ N, pii > 0. Stelling 1.27. Zij P irreducibel en aperiodiek en stel dat P een invariante verdeling π heeft. Zij λ een willekeurige distributie. Als (Xn )n≥0 Markov(λ, P ) is dan geldt P(Xn = j) → πj als n → ∞ voor alle j, en in het bijzonder, (n)

pij → πj als n → ∞ voor alle i, j. Bewijs. We zullen niet het hele bewijs in detail geven, maar een schets. Voor het gedetailleerde bewijs verwijzen wij naar het boek Markov Chains van Norris (1998) [5]. Het bewijs gaat in 3 stappen en maakt gebruikt van een koppelingsargument. Men definiëert een keten (Yn )≥0 die Markov(π, P ) is en onafhankelijk van (Xn )n≥0 . De eerste stap is te bewijzen dat voor een 13

vaste toestand b de Markov Chains elkaar in eindige tijd ontmoeten. Dus er wordt bewezen dat P(T = inf{n ≥ 1 : Xn = b = Yn } < ∞) = 1. Men doet dit door een Markovproces Wn = (Xn , Yn ) te maken met overgangskansen p˜(i,k)(j,l) = pij pkl en met beginverdeling µ(i,k) = λi πk . Daarna kan men door gebruik van de aperiodiciteit van P bewijzen dat P˜ irreducibel is en dat P˜ een invariante maat heeft. Met behulp stelling 1.25 is P˜ positief recurrent. Dit betekent met stelling 1.22 dat P(T < ∞) = 1. In de tweede stap definiëert men een Markovproces Zn = Xn als n < T en Zn = Yn als n ≥ T . Met behulp van de sterke Markoveigenschap voor (Wn )n≥0 op tijdstip T is (WT +n )n≥0 Markov(δ(b,b) , P˜ ) en onafhankelijk van (X0 , Y0 ), ..., (XT , YT ). Door ˜ T +n )n≥0 = (YT +n , XT +n )n≥0 Markov(δ(b,b) , P˜ ). gebruik van symmetrie is (W 0 0 0 We krijgen dan dat Wn = (Zn , Zn ) Markov(µ, P˜ ), waarbij Zn = Yn als n < T 0 en Zn = Xn als n ≥ T . In het bijzonder is dan (Zn )n≥0 Markov(λ, P ). In de derde stap kan men dan eenvoudig afleiden dat |P(Xn = j) − πj | ≤ 0 als (n) n → ∞, hetgeen betekent dat pij = πj . We hebben de nodige stellingen over Markov Chains gegeven en deels bewezen. In hoofdstuk 3 zullen er nog twee resultaten met betrekking tot Markov Chains volgen die nodig zijn voor Monte Carlo Markov Chain algoritmes.

14

Hoofdstuk 2 Bayesiaanse statistiek 2.1

Prior en posterior verdelingen

In dit hoofdstuk gaan we op de Bayesiaanse statistiek in. In de frequentistische statistiek wordt er aangenomen dat i.i.d. data X1 , ..., Xn gedetermineert worden door een vaste, maar onbekende parameter θ. Anders gezegd: Gegeven een model in de vorm van een verzameling kansmaten P = {Pθ : θ ∈ Θ}, zijn we op zoek naar de Pθ die de data zo goed mogelijk verklaart. Het doel is om dan met frequentistische inferentiële methoden een uitspraak te doen over de onbekende parameter θ. In de Bayesiaanse statistiek daarentegen wordt er geen onderscheid gemaakt tussen data X en parameters θ. Er wordt aangenomen dat θ een stochastische grootheid is op Θ. We willen dan een uitspraak doen over de verdeling op θ conditioneel op data X. Om dit te kunnen bewerkstellingen kiezen we een model voor de data conditioneel op θ. Het model is dan een verzameling kansmaten conditioneel op θ, dat wil zeggen P = {π(X|θ) : θ ∈ Θ}. Om uiteindelijk de conditionele verdeling op θ te vinden moeten we eerst kennis verzamelen over de verdeling van θ zelf. Dit doen we door een distributie θ 7→ π(θ) te geven. Deze verdeling noemen wij de prior. We kunnen met behulp van de regel van Bayes dan de posterior verdeling π(θ|X) afleiden. π(X|θ)π(θ) π(X, θ) =R π(X) π(X|θ) dΠ(θ) Θ π(X|θ)π(θ) =R . π(X|θ)π(θ) dθ Θ

π(θ|X) =

We zien dat we de posterior verdeling op θ gegeven data X op een constante na kennen, dit kan men ook schrijven als π(θ|X) ∝ π(X|θ) × π(θ). 15

Men kan ook π(X|θ) als de likelihoodfunctie bekijken. Deze noteren we dan met L(θ|X). We onthouden de volgende regel voor de posterior: π(θ|X) ∝ L(θ|X) × π(θ), posterior ∝ likelihood × prior . R In het algemeen zal het lastig zijn de integraal Θ π(X|θ)π(θ) dθ te berekenen als Θ bijvoorbeeld een hogere dimensionale ruimte is. We zullen in hoofdstuk 3 een uitweg vinden uit deze situatie. We kunnen dan zelfs zonder de normalisatieconstante te kennen direct uit de posterior trekken. Om wat gevoel voor een Bayesiaanse procedure te krijgen introduceren wij een voorbeeld van Rev. Thomas Bayes zelf (1763) [1],[4]. Voorbeeld 2.1. Gegeven een biljarttafel waarbij de langere zijde van de biljarttaffel lengte 1 heeft. Op deze tafel hebben wij n + 1 ballen, waarbij 1 bal wit is en de andere n zijn rood. Zij X ∈ [0, 1], de parameterruimte die de positie van de witte bal beschrijft. Laat (Y1 , ..., Yn ) de stochasten zijn die de positie tussen 0 en 1 van de rode ballen aangeeft. Het doel is, om op basis van K = aantal rode ballen dichter bij 0 dan de witte bal een aanpassing te geven van de a-priori uniform verdeelde positie van de witte bal.

Figuur 2.1: Een plaatje van de biljarttafel van Bayes. We zien een witte bal en n = 6 rode ballen. Er liggen K = 5 rode ballen dichter bij 0 dan de witte bal. Gegeven deze observatie, wat kunnen wij over de positie van de witte bal zeggen? Dan is (X, Y1 , ...Yn ) ∼ U [0, 1]n+1 verdeeld. De ballen zijn onafhankelijk op de tafel geplaatsd met als identieke distributies U [0, 1]. Wat weten we van X? Het enige wat we kunnen zeggen is dat de positie X = x tussen 0 en 1 16

ligt. Het beste wat we kunnen zeggen is dat X ∼ U [0, 1]. Deze nemen wij als prior. De observatie is P het aantal K rode ballen dichter bij 0 dan de witte bal. Dit betekent dat K = ni=1 1{Yi ≤ X}. Nadat wij deze observaties gemaakt hebben kunnen wij nog steeds niet zeggen wat P(X|K) is. Er geldt voor iedere i dat Yi en X onafhankelijk zijn en dus is P(Yi ≤ X|X = x) = P(Yi ≤ X) = P(Yi = x) = x. Dus kiezen we als model distributies P(Yi ≤ X|X = x) = x. Nu geldt voor iedere Yi dat deze boven de parameter x ligt of eronder. Dus voor iedere rode bal geldt als we willen vinden of deze dichter bij 0 ligt dan de witte bal dat we een Bernoulli experiment hebben met parameter x. P Omdat alle Yi i.i.d. zijn betekent dit dat K|X = x = ni=1 1{Yi ≤ X}|X = x binomiaal(x; n) verdeeld is. Dus P(K = k|X = x) = nk xk (1−x)n−k . Hieruit kunnen wij P(X = x|K = k) bepalen. P(X = x|K = k) = R = R1 0

P(K = k|X = x) P(K = k|X = x) dx x∈[0,1]

P(K = k|X = x) P(K = k|X = x) dx

= B(n, k)xk (1 − x)n−k ,

waarbij B(n, k) een normalisatie constante is. Deze normalisatieconstante moeten wij berekenen. We zien aan de vorm van X|K dat deze beta(k + 1, n − k + 1) verdeeld moet zijn. Dit betekent dat de normalisatieconstante B(n, k) = Γ(n + 2)/Γ(k + 1)Γ(k + 2) moet zijn. Wat betekent dit intu¨ıtief? Stel nou dat voor n = 1000 = K. Dat betekent dat alle rode ballen dichter bij 0 liggen dan de witte bal. Uit deze observatie mogen wij concluderen dat het zeer aannemelijk is dat de witte bal dicht bij 1 ligt. Omgekeerd als n = 1000 en K = 0 is het zeer waarschijnlijk dat de witte bal in de buurt bij 0 ligt. We hebben een plaatje gemaakt van de kansverdelingen voor n = 6 en bekijken de gevallen voor K ∈ {0, 1, ..., 6}.

17

7 K=0 K=1 K=2 K=3 K=4 K=5 K=6

6

5

4

3

2

1

0 0

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

1

Figuur 2.2: We zien de kansdichtheden van de posterior x|K bij behorende het voorbeeld van de biljarttafel van Bayes. Voor n = 6 laten we K = k lopen van 0 tot 6. Het moge met de intu¨ıtieve uitleg duidelijk zijn voor welke k wij welke kansdichtheid hebben. We sluiten deze sectie af met een definitie die belangrijk zal zijn voor hoofdstuk 3. Definitie 2.2. Zij (P, A ) een meetbaar model voor data X ∈ X . Zij M een verzameling verdelingen op (P, A ). De verzameling heet een geconjugeerde familie voor het model P, als de posterior gebaseerd op de prior uit de verzameling M weer een element van M is: Π ∈ M =⇒ Π(·|X = x) ∈ M voor alle x ∈ X . Als we een bepaalde keuze voor de prior hebben zodat de posterior ook in dezelfde familie van kansverdelingen ligt, is het mogelijk met hoge convergentiesnelheid in MATLAB de posterior te simuleren. In Hoofdstuk 3 zullen wij deze eigenschap gebruiken voor de Gibbs Sampler.

18

2.2

Inferenti¨ ele aspecten van de posterior

De posterior is het gereedschap dat een Bayesiaanse statisticus de mogelijkheid geeft de belangrijke onderdelen van de statistiek zoals schatting, het vinden van betrouwbaarheidsverzamelingen en toetsing te kunnen uitvoeren. We zullen dit aan hand van definities en voorbeelden duidelijker maken. Wij beginnen met een voorbeeld voor schatting. Definitie 2.3. Zij P een model met prior Π op Θ. Neem aan dat de posterior absoluut continu is ten opzichte van een σ eindige maat µ op P. We schrijven voor de µ dichtheid van Π(·, X) θ → π(θ|X). De maximum-a-posteriori ˆ (MAP) schatter θ(X) voor θ is gedefiniëerd als het punt in het model, waar de posterior verdeling zijn maximum aanneemt (gegeven dat er een maximum ˆ bestaat en uniek is). Dus θ(X) wordt zodanig gekozen dat: ˆ π(θ|X) = sup π(θ|X). θ∈Θ

Merk op dat de MAP-schatter een zwak punt heeft. Een andere keuze voor de dominiërende maat µ leidt tot een andere MAP-schatter. Voorbeeld 2.4. Stel wij hebben data X = (X1 , ..., Xn ) die onafhankelijk van elkaar π(X|θ) ∼ N (θ, σ 2 ) verdeeld zijn. Laat de prior voor θ π(θ) ∼ N (µ0 , σ0 2 ), waarbij µ0 , σ0 , σ bekend zijn. Dan is de posterior dichtheid π(θ|X) ∝ π(X|θ)π(θ). Omdat wij de posterior op een constante na kennen, is het voldoende om π(X|θ)π(θ) ten opzichte van θ te maximaliseren. π(X|θ)π(θ) (

−n 2

= (2πσ)

(2πσ0 )

− 21

1 exp − 2

2 n X Xi − θ σ

i=1

+

θ − µ0 σ0

2 !) .

Hieruit volgt dat wij 2 n X Xi − θ i=1

σ

+

θ − µ0 σ0

2

moeten maximaliseren. Door de afgeleide naar θ gelijk aan 0 te stellen verkrijgen wij ! n 2 X nσ 1 σ2 0 ˆ θ(X) = X + µ0 . i nσ0 2 + σ 2 n i=1 nσ0 2 + σ 2 19

Dit kan men ook schrijven als σ02 σ2 ˆ θ(X) = X n + 2 µ0 . σ02 + n−1 σ 2 nσ0 De MLE (maximum likelihood estimator) schatter voor θ is X n . We zien dat ˆ er sprake is van een n−1 -proportionele bias term die de MAP-schatter θ(X) richting µ0 drijft. Deze bias verdwijnt echter asymptotisch: σ02 σ2 lim X n + 2 µ0 = X n . n→∞ σ02 + n−1 σ 2 nσ0 en we zien dat in dit voorbeeld de MAP schatter convergeert naar de MLE schatter. We merken op dat er ook andere methoden zijn voor schatting zoals de posterior mean of de posterior median. Voor verdere voorbeelden zie de syllabus voor Bayesian statistics (2011) [4] en het boek the Bayesian Choice (2007) [1]. We zullen nu een manier aangeven hoe een Bayesiaan betrouwbaarheidsverzamelingen met behulp van de posterior kan vinden. Definitie 2.5. Zij (Θ, G ) een meetbare ruimte die een model Θ → P : θ → Pθ voor data X parametrizeert, met een prior Π : G → [0, 1]. Kies een α ∈ (0, 1). Zij D ∈ G een deelverzameling van Θ. Dan noemen we D een level-α credible set, als Π(ϑ ∈ D|X) ≥ 1 − α. Definitie 2.6. We definiëren voor iedere k ≥ 0 de level-set: D(k) = {θ ∈ Θ : π(θ|X) ≥ k}. Definitie 2.7. Zij (Θ, G ) een meetbare ruimte die een model Θ → P : θ → Pθ voor data X ∈ X parametrizeert met een prior Π : G → [0, 1]. Neem aan dat de posterior gedominiëerd wordt door een σ-eindige maat µ op (Θ, G ), met dichtheid π(·|X). Kies een α ∈ (0, 1). Een level-α HPD-credible set (HPD=highest posterior density) voor ϑ is de deelverzameling Dα = D(kα ), waarbij, kα = sup{k ≥ 0 : Π(ϑ ∈ D(k)|X) ≥ 1 − α}. Voorbeeld 2.8. Zij X = (X1 , ..., Xn ) een i.i.d. steekproef uit N (θ, σ 2 ) met θ ∈ R onbekend en σ 2 bekend. Laat de prior π(θ) ∼ N (µ0 , σ02 ), waarbij µ0 en σ02 bekend zijn. De posterior is dan π(θ|X) ∼ N (µp , σp2 ), waarbij 20

2

σ2 σ2

nσ 2

2 0 0 µp = nσ2σ+σ2 µ0 + nσ2 +σ 2 X n en σp = nσ 2 +σ 2 . Een level-α HPD credible set 0 0 0 wordt gegeven door n o √ 2 D = θ ∈ Θ : e−(µp −θ)/2σp ≥ kα 2πσp o n √ √ 1 2 = θ ∈ Θ : |θ − µp | ≤ 2σp [− log(kα 2πσp )] .

We willen nu dat π(θ|X)(D) = 1 − α. Met een beetje rekenwerk komt dit erop neer dat √ √ 1 α , 2σp [− log(kα 2πσp )] 2 = σp Φ−1 1 − 2 hetgeen betekent dat h α α i , µp + σp Φ−1 1 − D = µp − σp Φ−1 1 − 2 2 " # σσ0 σσ α α 0 = µp − p 2 , µp + p 2 Φ−1 1 − Φ−1 1 − 2 2 nσ0 + σ 2 nσ0 + σ 2 een HPD-credible set voor θ is. Hypothese-toetsing in Bayesiaanse context lijkt veel op classificatie lijkt en dit is aanzienlijk eenvoudiger dan het populaire Neyman-Pearson-paradigma uit de frequentistische statistiek. Ten slotte geven wij een voorbeeld van toetsing in Bayesiaanse context. Definitie 2.9. Zij (Θ, G ) een meetbare ruimte die een model Θ → P : θ → Pθ voor data X ∈ X parametrizeert met een prior Π : G → [0, 1]. Zij {Θ0 , Θ1 } een partitie van Θ zdd. Π(Θ0 ) > 0 en Π(Θ1 ) > 0. We schrijven voor de posterior odds: Π(Θ0 |X) , Π(Θ1 |X) en voor de prior odds Π(Θ0 ) . Π(Θ1 ) De Bayesfactor met voorkeur voor de deelverzameling Θ0 is B=

Π(Θ0 |X) Π(Θ1 ) . Π(Θ1 |X) Π(Θ0 ) 21

Dit betekent voor een (subjectieve) Bayesiaan dat als in onze definidat we de voorkeur geven aan de verzameling Θ0 en als

0 |X) tie Π(Θ > 1, Π(Θ1 |X) Π(Θ0 |X) < 1, dat Π(Θ1 |X)

we de voorkeur geven aan Θ1 . Dit zullen we aan hand van een voorbeeld duidelijker maken: Voorbeeld 2.10. We nemen als voorbeeld de posterior verdeling van het biljardprobleem van Bayes. De posterior heeft dichtheid π(x|K = k) = Γ(n + 2)/Γ(k + 1)Γ(k + 2)xk (1 − x)n−k . Laat het aantal rode ballen n = 10 zijn en laat k = 8 zijn. We kiezen Θ0 = [0, 21 ] en Θ1 = ( 12 , 1]. Dit komt er op neer om H0 : x ≤ 1/2 tegen H1 : x > 1/2 te toetsen. De Bayesfactor met voorkeur voor Θ0 is R1 8 2 x (1 − x)2 dx 1/2 67/1013760 67 Π(Θ0 |K) Π(Θ1 ) 0 = R1 = = . B= 8 2 Π(Θ1 |K) Π(Θ0 ) 1/2 1981/1013760 1981 1 x (1 − x) dx 2

67 Omdat 1981 < 1, concluderen wij dat het aannemelijker is dat x ∈ Θ1 . We verwerpen H0 en accepteren H1 .

22

Hoofdstuk 3 Monte Carlo Markov Chains Wij hebben nu het gereedschap verzameld om algoritmes te introduceren die wij Monte Carlo Markov Chain algoritmes noemen. (Deze worden vaak ook wel afgekort met MCMC.) We construeren hierbij Markov Chains die als limietverdeling onze doeldistributie zullen krijgen. We merken hierbij op dat wij meestal met toestandsruimtes werken die overaftelbaar zijn. Dit heeft tot gevolg dat wij alleen een discrete benadering van ons probleem kunnen geven. Wij zullen in dit hoofdstuk twee algoritmes geven en aan het eind van ons hoofdstuk de algoritmes toepassen voor bekende dichtheden zoals de Cauchyverdeling, maar ook voor posterior verdelingen.

3.1

Symmetrisch Metropolis Hastings algoritme

We geven eerst een definitie en een stelling. Definitie 3.1. Een stochastische matrix P en een maat π zijn in detailed balance, als πi pij = πj pji ∀i, j ∈ I. Stelling 3.2. Als P en π in detailed balance zijn, dan is π een invariant voor P . Bewijs. (πP )i =

X

πj pji =

j∈I

X j∈I

23

πi pij = πi .

Voor we het algoritme geven, geven wij uitleg over notatie. De toestandsruimte van de Markov Chain die wij construeren is gelijk aan de parameterruimte Θ. Verder nemen wij aan dat de overgangskansen van de proposaldistributie symmetrisch zijn, dat q(x|y) = q(y|x). We construeren een Markov Chain als volgt. 1. Kies een startwaarde θ0 . 2. Gegeven θn , trek vervolgens een proposal Yn+1 ∼ q(θi |θj ) = q(θj |θi ). o n n+1 ) . 3. Bereken de acceptatieratio α(θn , Yn+1 ) = min 1, π(Y π(θn ) 4. Zij U ∼ U [0, 1]. 5. Accepteer proposal Yn+1 als U < α. 6. Als we accepteren zetten we θn+1 = Yn+1 en θn+1 = θn elders. 7. Herhaal stap 2 t/m 6! We merken op voor het geval dat π(θn ) = 0 wij α(θn , Yn+1 ) = 1 zetten. De keuze voor de acceptatiekans α komt neer op de volgende propositie: Propositie 3.3. Het symmetrische Metropolis Hastings algoritme genereert een Markov Chain (θn )n≥0 die invariant is ten opzichte van π(·). Bewijs. Wij moeten de detailed balance vergelijking nagaan: We moeten voor alle θi , θj ∈ Θ laten zien of π(θi ) q(θj |θi ) α(θi , θj ) = π(θj ) q(θi |θj ) α(θj , θi ) Nu is q(θi |θj ) = q(θj |θi ), dit heeft tot gevolg dat π(θi ) q(θi |θj ) α(θi , θj ) π(θj ) = π(θi ) q(θj |θi ) min 1, π(θi ) = q(θj |θi ) min {π(θi ), π(θj )} π(θi ) = q(θj |θi ) π(θj ) min ,1 π(θj ) = π(θj ) q(θj |θi ) α(θj , θi ). Het symmetrische Metropolis Hastings algoritme voldoet aan de detailed balance vergelijking. Dus is π invariant voor de keten. 24

Opmerking 3.4. Als wij uit een posterior met dichtheid π(θ|X) willen trekken zijn er drie ingrediënten nodig. Wij hebben de dichtheid van de prior π(θ) nodig en we hebben een model nodig voor onze data. In onze voorbeelden hebben wij deze ingediënten gegeven. Verder moeten wij een symmetrische propsaldistributie voor het algoritme kiezen. We zullen daarvoor de normale verdeling gebruiken. Andere symmetrische distributies zijn natuurlijk ook mogelijk. Herinner dat de dichtheid van de posterior geschreven kan worden als π(X|θ)π(θ) R . π(X|θ)π(θ) dθ Θ R Omdat wij Θ π(X|θ)π(θ) dθ niet kennen zouden wij voor deze integraal Monte Carlo integratie kunnen toepassen. Dit hoeft echter niet: Wij kunnen namelijk meteen uit de posterior simuleren, omdat voor alle θi , θj ∈ Θ geldt dat R π(X|θ)π(θ) dθ π(X|θ)(θj )π(θj ) π(θ|X)(θj ) =R · Θ π(θ|X)(θi ) π(X|θ)π(θ) dθ π(X|θ)(θi )π(θi ) Θ π(X|θ)(θj )π(θj ) = , π(X|θ)(θi )π(θi ) en zowel de prior als ook het model zijn in θi respectievelijk θj bekend. Hieruit volgt voor posteriors dat de acceptatiekans π(X|θ)(θj )π(θj ) α(θi , θj ) = min 1, π(X|θ)(θi )π(θi ) is. We kunnen dus met behulp van het symmetrische Metropolis Hastings algoritme direct uit de posterior trekkingen doen. We hebben nu een algoritme gezien dat het mogelijk maakt een posterior verdeling op basis van een keuze voor het model en een prior te benaderen. Een nadeel van het symmetrische Metropolis-Hastings algoritme is dat een getrokken proposalwaarde een bepaalde acceptatiekans heeft. Het is mogelijk dat waarden gegenereerd door de proposalverdeling niet geaccepteerd worden. We zullen daarom een ander algoritme bekijken dat helemaal geen acceptatievoorwaarde heeft. We introduceren in de volgende sectie de Gibbs Sampler.

3.2

Gibbs Sampler

Voor dat wij het algoritme van de Gibbs sampler geven, zijn er nog belangrijke notatiekwesties te behandelen. Zij X ⊂ Rm en zij x = (x1 , ..., xm )T 25

een vector in X . Zij {π(x) : x ∈ X } een verdeling op X . We noteren de stochastische vector X = (X1 , ..., Xm )T zo dat P(X = x) = P(X1 , = x1 , ..., Xm = xm ) = π(x1 , ..., xm ), voor alle x ∈ X . Voor de marginale dichtheden van Xi schrijven we πXi zo dat P(Xi = xi ) = πXi . Voor een deelverzameling T ⊆ M = {1, 2, ..., m} is XT de vector met de componenten Xi , waarvoor geldt dat i ∈ T en X−T is de vector met de componenten Xi , i 6∈ T . De marginale kans van XT schrijven wij πXT en voor de marginale kans van X−T schrijven wij πX−T . Met deze notatie kunnen wij nu handig conditionele kansen weergeven. Het is duidelijk dat P(XT = xT |X−T = x−T ) = πXT |X−T =x−T (xT ). Zij T een vast gekozen deelverzameling van M. We construeren nu een Markov Chain (Xn ) afhankelijk van T als volgt. Zij X de toestandsruimte 0 0 van de Markov Chain. Gegeven X = x zij X 0 = (X1 , ..., Xm ) de volgende 0 0 toestand in de keten. Dan is gegeven X = x met X−T = x−T en wordt XT 0 onafhankelijk getrokken met kans πXT |X−T =x−T (xT ). Dan vormt (Xn ) een Markov Chain met overgangskansen ( 0 0 als x−T 6= x−T pxx0 = 0 0 πXT |X−T =x−T (xT ) als x−T = x−T 0

voor alle x, x ∈ X . 0

Stelling 3.5. Als X de distributie π heeft dan heeft X ook distributie π . Bewijs. We ontlenen het bewijs uit ”Monte Carlo Markov Chain Simulation” (2011) van Bert van Es [3]. Er geldt: X X 0 0 0 0 P(X ∈ C) = P(X = x) = P(XT = xT , X−T = x−T ) x∈C

=

X

=

X

x∈C 0

0

0

P(XT = xT |X−T = x−T ) · P(X−T = x−T )

x∈C

πX 0 |X 0 T

−T =x−T

(xT ) · πX 0 (x−T ) −T

x∈C

=

X

πXT |X−T =x−T (xT ) · πX−T (x−T ) =

x∈C

X x∈C

Dus π is een invariante distributie voor de keten. 26

πX (x) = π(C)

Opmerking 3.6. Het feit dat de keten een invariante verdeling heeft betekent niet dat de keten ook daadwerkelijk naar de invariante verdeling convergeert. De eisen voor convergentie zijn namelijk irreducibiliteit en aperiodiciteit. De aperiodiciteit volgt uit het feit dat als Xn = xn en we willen een stap doen in de Markov Chain, dat pxn xn > 0. Stel maar dat pxn xn = 0, dan hadden we xn bij de vorige stap niet kunnen bereiken. Het probleem wat zich voor doet is dat de keten niet irreducibel hoeft te zijn. Als namelijk de keten een willekeurige begintoestand heeft, dan is de kansverdeling van de −T -componenten hetzelfde als de kansverdeling onder de startverdelingen. Voor een manier om de Gibbs Sampler toch irreducibel te krijgen kiezen we voor iedere stap die we in de keten willen maken een andere deelverzameling T . Om hierover meer te weten te komen verwijzen wij naar de literatuur [3]. Onze keuze voor T zal T = {i} zijn. De zo verkregen keten wordt ook wel single site Gibbs sampler genoemd. De overgangskansen van de keten worden dan gegeven door: ( 0 0 als x−i 6= x−i pxx0 = 0 0 πXi |X−i =x−i (xi ) als x−i = x−i De vraag is hoe de Gibbs Sampler er voor posterior verdelingen uit ziet. Om dit te beantwoorden zullen we de single site Gibbs sampler gebruiken: Zij θ = (θ1 , ..., θm )T de parameter vector, L(X|θ) de likelihood en π(θ) de prior voor θ. We hebben in Hoofdstuk 2 gezien dat π(θi |θj , i 6= j, X) ∝ L(X|θ)π(θ). We geven het algoritme van de single site Gibbs sampler voor posteriors: (0)

(0)

1. Zij t = 0 en kies een startwaarde θ(0) = (θ1 , ..., θm ) en een stoptijd T > 0. 2. Genereer de componenten van θ als volgt: (t+1)

uit π(θ1 |θ2 , ..., θm , X)

(t+1)

uit π(θ2 |θ1

(t+1)

uit π(θm |θ1

• trek θ1

• trek θ2

(t)

(t+1)

(t)

(t)

(t)

, θ3 , ..., θm , X)

• ... • trek θm

(t+1)

(t+1)

, ..., θm−1 , X)

3. Zet t := t + 1. Als t < T herhaal de vorige stap. We zien dat wij de verdeling van de parameter θ gegeven X kunnen simuleren door uit de conditionele dichtheden op π(θi |...) te simuleren. Dit is alleen nuttig als we weten hoe de posterior eruit ziet. In sommige gevallen is het niet mogelijk de conditionele verdelingen te vinden en kunnen wij beter het symmetrische Metropolis Hastings algoritme gebruiken. We hebben 27

in Hoofdstuk 2 echter een definitie van een geconjugeerde familie gegeven. Dit betekent dat als gegeven een model en een prior de posterior in dezelfde familie van kansverdelingen zit als de prior, dat wij de Gibbs Sampler de voorkeur zullen geven. Dit is omdat de Gibbs Sampler geen acceptatievoorwaarde heeft. We geven in de volgende sectie voorbeelden voor het uitvoeren van het symmetrische Metropolis Hastings algoritme en de Gibbs Sampler. We gebruiken de MCMC algoritmes daarbij ook voor bekende kansdichtheden zoals bijvoorbeeld de Cauchyverdeling maar zullen ook voorbeelden geven van onbekende posterior verdelingen.

3.3

Voorbeelden

In deze sectie werken we voorbeelden voor het symmetrisch Metopolis Hastings algoritme en de Gibbs Sampler uit. Voorbeeld 3.7. Stel wij willen trekkingen doen uit de Cauchyverdeling met 1 dichtheid f (x) = π1 · 1+x 2 . Dit kan met het symmetrische Metropolis Hastings algoritme. We kiezen als proposaldichtheid de normale verdeling. We testen het algoritme voor verschillende varianties σ van de proposalverdeling. We geven het algoritme.

1. Zij n = 0 en kies een willekeurige startwaarde Xn . 2. Zij Y de proposal q(Y |Xn ) ∼ N (Xn , σ) verdeeld 3. Bereken de ratio α(Xn , Y ) f (Y )q(Xn |Y ) ,1 = min f (Xn )q(Y |Xn ) 1 + Xn 2 = min ,1 1+Y2

f u n c t i o n [X] = SymMetropolisHastingsCauchy (N, sigma , x0 ) X( 1 )=randn ; f o r i =2:N Y=normrnd (X( i −1) , sigma ) ; % p r o p o s a l voor Y a l p h a=min ( (1+(X( i −1)−x0 ) ˆ 2 ) /(1+( Y−x0 ) ˆ 2 ) , 1 ) ; %r a t i o u=rand ; i f u
4. Zij U ∼ U [0, 1]. Accepteer als U < α, verwerp elders en zij n := n + 1, herhaal stap 2 t/m 4.

28

4000

700

3500

600

3000

500

2500 400 2000 300 1500 200

1000

100

500 0 −30

−20

−10

0

10

20

0

30

0

0.2

0.4

0.6

0.8

1

1.2

1.4

Figuur 3.1: Het histogram links is verkregen met behulp van het symmetrische Metopolis Hastings algoritme met als proposaldichtheid een normale verdeling met variantie σ = 2. Rechts zien we een histogram dat door dezelfde methode verkregen wordt, maar dan met σ = 0.01. De σ in ons algoritme geeft aan met welke stapgroote wij een nieuwe proposal trekken. Als we bijvoorbeeld een hele kleine σ van 0.01 kiezen, zullen we maar in een klein gebied rond onze huidige toestand θn een nieuwe proposal trekken. Een voordeel kan zijn dat een proposal in de buurt van θn een hoge acceptatiekans heeft, maar een nadeel is wel dat wij door hele kleine stappen te maken niet de hele toestandsruimte te zien krijgen. We zien dan voor het geval σ = 0.01 dat de Cauchyverdeling nog niet goed benaderd wordt. Kiezen we dan σ = 2 kan de acceptatiekans van nieuwe proposalwaarden kleiner zijn, maar heeft deze keuze van σ wel als voordeel dat we sneller een groter gedeelte van de toestandsruimte gezien hebben. Dit is goed te zien in figuur 3.2.

29

2

4.6

0

4.4

−2 4.2

−6

(θn)n≥ 0

(θn)n ≥ 0

−4

−8

4

3.8

−10 3.6 −12 3.4

−14 −16 500

510

520

530

540

550 560 iteraties n

570

580

590

3.2 500

600

510

520

530

540

550 560 iteraties n

570

580

590

600

Figuur 3.2: We zien links het pad van de Markov Chain voor σ = 2 en rechts het pad van de Markov Chain voor σ = 0.01. Als we deze twee vergelijken, valt op dat de Markov Chain links grotere sprongen maakt door de parameterruimte, maar daarvoor vaker blijft staan. De Markov Chain rechts beperkt zich tot een veel kleiner gedeelte van de parameterruimte, maar accepteert vaker nieuwe proposals. We zullen nu het symmetrische Metropolis Hastings algoritme gebruiken voor de benadering van een posterior. Hierbij zullen we gebruiken dat wij uit de Cauchyverdeling met willekeurige locatieparameter kunnen trekken. Voorbeeld 3.8. We bekijken het voorbeeld 6.1.1 uit het boek ”The Bayesian Choice” [1]. Gegeven een steekproef X1 , ..., Xn uit de Cauchyverdeling C(θ, 1) met locatieparameter θ en gegeven een prior op θ die N (µ, σ 2 ) verdeeld is, waarbij µ en σ bekend zijn. We krijgen voor de posterior verdeling dat: 2 /2σ 2

π(θ|X1 , ..., Xn ) ∝ e(θ−µ)

n Y [1 + (Xi − θ)2 ]−1 i=1

Omdat men de posterior op een constante na kent is het handig het symmetrische Metropolis Hastings algoritme te gebruiken. Het vorige voorbeeld maakt het mogelijk uit een Cauchyverdeling met willekeurige locatieparameter trekken. We voeren het symmetrische Metopolis Hastings algoritme uit.

30

f u n c t i o n [ t h e t a ] = s y m m e t r o p o l i s h a s t i n g s p o s t e r i o r (mu, sigma , N1 , N2) t h e t a ( 1 ) =0; produkt =1; f o r i =2:N2 X= SymMetropolisHastingsCauchy (N1 , 0 . 2 , t h e t a ( i −1) ) ; % data Y=randn ; % p r o p o s a l v o o r Y f o r j =1:N1 produkt=produkt ∗( (1 +(X( j )−t h e t a ( i −1) ) ˆ 2 ) /(1+(X( j )−Y) ˆ 2 ) ) ; end a l p h a=min ( exp (( −(Y−mu) ˆ2+( t h e t a ( i −1)−mu) ˆ 2 ) / ( 2 ∗ sigma ˆ 2 ) ) ∗ produkt , 1 ) ; u=rand ; i f u
Figuur 3.3: Het symmetrische Metopolis Hastings algoritme in een MATLAB code voor de posterior

3500

3000

1800 1600

3000

2500 1400

2500 2000

1200

2000

1000 1500 800

1500 1000

600

1000 400 500

500

0 −6

200

−4

−2

0

2

4

6

0 −3

−2

−1

0

1

2

3

0 −0.3

−0.2

−0.1

0

0.1

0.2

0.3

0.4

0.5

0.6

Figuur 3.4: Histogram van simulatie met N = 40000 trekkingen gegenereerd door het symmetrische Metropolis Hastings algoritme met als proposalverdeling de normale verdeling. Links hebben we voor de variantie van de proposal σk = 4 gekozen. In het midden hebben we σk = 0.4 gekozen. Rechts hebben we σk = 0.04 gekozen. De beste benadering van de posterior geeft het histogram in het midden. Ook hier zien we weer dat verschillende waarden van σk tot snellere of langzamere convergentie tot de posterior kunnen leiden. De reden hiervoor is dezelfde als bij het voorbeeld van de Cauchy verdeling 3.7. We hebben voorbeelden voor het symmetrische Metropolis Hastings algoritme gegeven. Het symmetrische Metropolis Hastings is handig om posterior verdelingen die niet in gesloten vorm beschikbaar zijn discreet te benaderen. We geven nu voorbeelden van kansverdelingen waar voor wij de Gibbs Sampler zullen gebruiken. 31

Voorbeeld 3.9. Stel wij willen steekproeven genereren uit de verdeling f (θ1 , θ2 ) die N (µ, Σ) verdeeld is, met µ = (1, 2)T en ! 1 ρ Σ= ρ 1 We weten voor de conditionele kansdichtheden dat f (θ1 |θ2 ) ∼ N (µ1 + ρ(θ2 − µ2 )), 1 − ρ2 ) en f (θ2 |θ1 ) ∼ N (µ2 + ρ(θ1 − µ1 )), 1 − ρ2 ). We voeren de Gibbs Sampler uit:

1. Kies t = 0 en T en zij σ =

p 1 − ρ2 .

2. Kies een startwaarde voor (θ1 (t+1)

3. trek θ1

(t)

θ2 ) = (0, 0)T . (t+1)

(t)

uit f (θ1 |θ2 ) en trek θ2

(t)

(t+1)

uit f (θ2 |θ1

)

4. Herhaal tot t = T .

f u n c t i o n [ t h e t a ] = Gibbs (T) Y = [ 1 ; 2 ] ; %mu = [ 1 ; 2 ] rho = 0 . 0 1 ; %k i e s een g e t a l rho ; Hiermee maken we s t a p p e n i n de markovketen sigma = s q r t ( 1 − rho ˆ 2 ) ; %d u i d e l i j k ! t h e t a ( 1 , : ) = [ 0 0 ] ; %s t a r t w a a r d e f o r t = 2 :T mu1 = Y( 1 ) + rho ∗ ( t h e t a ( t −1 ,2) − Y( 2 ) ) ; %s t a p p e n v o o r de eerste coordinaat t h e t a ( t , 1 ) = normrnd (mu1 , sigma ) ; % sample u i t n o r m a l e v e r d e l i n g met p a r a m e t e r s mu1 en sigma mu2 = Y( 2 ) + rho ∗ ( t h e t a ( t , 1 ) − Y( 1 ) ) ; %s t a p p e n v o o r de tweede coordinaat t h e t a ( t , 2 ) = normrnd (mu2 , sigma ) ; %sample u i t normale v e r d e l i n g met p a r a m e t e r s mu2 en sigma end ; theta contour ( theta ( : , 1 ) , theta ( : , 2 ) ) end

We voeren de simulatie in MATLAB uit en krijgen de volgende resultaten.

32

7

5

6 4 5 4

3

θ2

θ2

3 2

2 1

1

0 0 −1 −2 −3

−2

−1

0

1

θ1

2

3

4

5

−1 −2

6

−1

0

1 θ1

2

3

4

Figuur 3.5: We zien links een MATLAB plot van de vector (θ1 , θ2 ) na T = 20000 iteraties. Rechts zien we het pad dat de Markov Chain aflegt na de eerste 300 iteraties. We zien hoewel wij in (0, 0) begonnen zijn, dat de getrokken punten zich centreren rond (θ1 , θ2 ) = (1, 2).

6 0.15 4 0.1 2 0.05 0

−2

−1

−2 0

1

2

3

4

θ2

5

θ1

Figuur 3.6: Een 3d-plot van de multivariaat normale verdeling uit voorbeeld 3.9 met 300 datapunten (θ1 , θ2 ) verkregen door de Gibbs Sampler.

33

6 5 4

θ2

3 2 1 0 −1 −2 −2

−1

0

1

θ1

2

3

4

5

Figuur 3.7: We zien een contourplot verkregen door de Gibbs Sampler uit voorbeeld 3.9 voor 20000 punten. Ook hier zien we duidelijk dat de punten die we met de Gibbs Sampler trekken zich rond (θ1 , θ2 ) = (1, 2) centreren. Voorbeeld 3.10. Wij geven nu een voorbeeld voor een Gibbs Sampler die trekking uit een posterior doet. Gegeven een model X|p ∼ bin(n, p). Zij de π(p) de prior voor p met dichtheid Beta(α, β), waarbij α en β bekend zijn. De dichtheid is dan π(p|X) ∝ pX (1 − p)n−X pα−1 (1 − p)β−1 = pα+X−1 (1 − p)n−X+β−1 . We concluderen dat de posterior in de famile van Betaverdelingen ligt. Dit is een voorbeeld van een geconjugeerde familie zoals wij dat in hoofdstuk 2 gedefiniëerd hebben. De posterior is Π(·, X) = Beta(α +X −1, n−X +β −1) verdeeld. Omdat MATLAB een algoritme heeft om uit de beta verdeling te trekken is mogelijk met de Gibbs Sampler meteen uit de posterior te trekken.

34

f u n c t i o n [ t h e t a ] = G i b b s c o n j u g a t e 1 (T, alpha , beta , n ) p=0.5; t h e t a ( 1 )=p ; t h e t a a c c e n t ( 1 ) =0; f o r t = 2 :T t h e t a a c c e n t ( t )=b i n o r n d ( n , t h e t a ( t −1) ) ; t h e t a ( t )=b e t a r n d ( a l p h a+t h e t a a c c e n t ( t −1) , b e t a +(n− t h e t a a c c e n t ( t −1) ) ) ; end ; theta h i s t ( theta , 1 5 5 ) end

4000 3500 3000 2500 2000 1500 1000 500 0

0

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

1

Figuur 3.8: Simulatie van de posterior die in dezelfde familie van kansverdelingen als de prior ligt. We verkrijgen een histogram door de Gibbs Sampler uit te voeren voor T = 100000, α = 2, β = 5, n = 2. Voorbeeld 3.11. Dit voorbeeld is uit het boek ”The Bayesian Choice” van Christian Robert [1]. Bekijk (λ, µ) ∈ N ∪ {0} × [0, 1] en de posterior n λ+α−1 π(λ, µ|X) ∝ µ (1 − µ)n−λ+β−1 , λ waarbij α en β dataafhankelijke variabelen zijn. De conditionele dichtheden die nodig zijn voor de implementatie van de Gibbs Sampler zijn λ|X, µ ∼ bin(n, µ) en µ|X, λ ∼ Beta(α + λ, β + n − λ). We kunnen met behulp van de Gibbs Sampler nu π(λ, µ|X) simuleren.

35

1. Kies θ(0) = (λ(0) , µ(0) ) = (0, 0). Kies n, α, β vast. 2. Zij t = 1. Trek λ(t) uit π(λ|X, µ(t−1) ) en trek daarna µ(t) uit π(µ|X, λ(t) ). 3. Herhaal tot t = T .

1

1

0.9

0.9

0.8

0.8

0.7

0.7

0.6

0.6

0.5

0.5

µ

µ

f u n c t i o n [ t h e t a ] = Gibbs4 (T, alpha , beta , n , waarnemingbegin , waarnemingeinde ) lambda =0; mu=0; t h e t a ( 1 , : ) =[ lambda mu ] ; %b e g i n w a a r d e s f o r t = 2 :T lambda = b i n o r n d ( n , mu) ; t h e t a ( t , 1 )=lambda ; mu = b e t a r n d ( a l p h a+t h e t a ( t , 1 ) , b e t a+n−t h e t a ( t , 1 ) ) ; t h e t a ( t , 2 )=mu ; % end p l o t ( t h e t a ( waarnemingbegin : waarnemingeinde , 1 ) , t h e t a ( waarnemingbegin : waarnemingeinde , 2 ) , ’ . ’ ) end

0.4

0.4

0.3

0.3

0.2

0.2

0.1 0

0.1

0

1

2

3

4

5 λ

6

7

8

9

0

10

0

10

20

30

40

50 λ

60

70

80

90

100

Figuur 3.9: Links zien we een plot van (λ, µ) verkregen door de Gibbs Sampler. We hebben de plot verkregen voor de parameters T = 20000, α = β = 0.1, n = 10. Rechts hebben we een andere plot van (λ, µ) verkregen door de Gibbs Sampler met parameters T = 20000, α = 0.5, β = 0.5, n = 100.

36

1 0.9

0.8

0.8

0.7

0.7

0.6

0.6

0.5

0.5

µ

µ

1 0.9

0.4

0.4

0.3

0.3

0.2

0.2

0.1 0

0.1

0

1

2

3

4

5 λ

6

7

8

9

0

10

0

10

20

30

40

50 λ

60

70

80

90

100

Figuur 3.10: Links hebben wij een pad gemaakt van de Markov Chain die gegenereerd wordt door de Gibbs Sampler uit het vorige figuur links. Rechts zien we een soortgelijk pad maar dan voor de Gibbs Sampler uit het vorige figuur rechts.

37

Hoofdstuk 4 Conclusie en discussie 4.1

Samenvattend

De vraag die we ons stelden in het begin van deze scriptie was of het mogelijk was de posterior verdeling te simuleren. Om deze vraag te beantwoorden moesten wij enige kennis verzamelen. Zo zijn we in hoofdstuk 1 begonnen met de basis. We hebben ons geconcentreerd op kansverdelingen wiens verdelingsfunctie een strict stijgende continue functie is. Daar hebben we het inversiealgoritme op toegepast. Het maakte ons mogelijk de exponentiële verdeling te simuleren. Vanaf daar bekeken wij de verwerpingsmethode. We gebruikten dat we met behulp van de inversiemethode de exponentiële verdeling konden simuleren om in combinatie met de verwerpingsmethode uit een normale verdeling met willekeurige parameters te kunnen trekken. Daarna hebben we ons in de theorie van de Markov Chains verdiept. We hebben belangrijke stellingen af kunnen leiden om uiteindelijk naar het einddoel in hoofdstuk 1 te komen. We lieten zien dat als een Markov Chain aperiodiek en irreducibel is, dan convergeert de chain naar zijn invariante maat. Deze stelling zou ons het idee geven om een Markov Chain zo te construeren, dat deze een invariante maat heeft die gelijk is aan de te simuleren verdeling. Omdat wij ge¨ınteresseerd zijn in het simuleren van posterior verdelingen, hebben wij in hoofdstuk 2 de Bayesiaanse statistiek ge¨ıntroduceert. Allereerst moesten wij begrijpen wat een posterior verdeling is en hebben we de problematiek van de posterior bepaald. De posterior is namelijk alleen op normalisatie na beschikbaar. Door naar de inferentiële aspecten van de Bayesiaanse statistiek te bekijken, hebben we inzicht kunnen krijgen in de waarde van de posterior voor een Bayesiaanse statisticus. Dit motiveerde ons om de 38

Monte Carlo Markov Chains te gebruiken om de posterior te benaderen. In hoofdstuk 3 hebben wij twee MCMC algoritmes, het symmetrische Metopolis Hastings algoritme en de Gibbs Sampler gegeven voor willekeurige verdelingen. Daarna hebben wij ons geconcentreerd op het toepassen van deze algoritmes op een willekeurige posterior verdeling. Uiteindelijk waren we in staat de praktische invulling met behulp van een computer en MATLAB te geven. We verzonnen ten dele eigen voorbeelden en bekeken voorbeelden uit The Bayesian Choice [1]. Voor al deze voorbeelden hebben we een MATLAB implementatie van de MCMC algoritmes gemaakt zo dat wij als invariante verdeling van de keten onze doeldistributie hadden. De resultaten waren daarbij meestal goed. Alhoewel wij ons er niet van konden verzekeren dat onze ketens irreducibel en aperiodiek waren, leek het erop dat onze ketens naar de gewenste verdeling convergeerden. Een klein probleem vonden we bij het symmetrische Metropolis Hastings algoritme. Een verandering van de parameterwaarde van de bekende proposalverdeling kon voor snellere of langzamere convergentie naar de doeldistributie zorgen. Een voorbeeld is de posterior uit voorbeeld 3.8. Op basis van ervaringen door bijvoorbeeld het algoritme 100000 keer te herhalen en te experimenten kon men een indruk krijgen hoe de invariante maat eruit moest zien. Omdat 100000 herhalingen een tijdrovende bezigheid is, wilden we kijken of we een dergelijk resultaat konden krijgen door minder vaak het algoritme te herhalen. Dit deden we door de variantieparameter van de proposalverdeling op verschillende manieren te kiezen. Zo was het mogelijk aan de ene kant mogelijk dat we voor te grote waarden van de varianties men na 40000 keer te herhalen een vervalst resultaat verkregen. Dit impliceerde dat de verdeling nog niet naar de invariante maat convergeerde. Men kon uiteindelijk door te experimenteren het resultaat optimalizeren, zo dat we na 40000 herhalingen een resultaat hadden wat sneller naar de invariante maat convergeerde. Ten slotte hebben ons de Monte Carlo Markov Chains de mogelijkheid gegeven om problemen in de Bayesiaanse statistiek, het vinden van de posterior, numeriek te benaderen. Met behulp van MCMC kunnen we zelfs bepaalde verwachtingen, varianties en hogere momenten, door Monte Carlo integratie te gebruiken, schatten. Dit was voorheen niet mogelijk.

39

4.2

Resterende vragen

Er resteren wel nog enkele vragen. We kunnen Markov Chains construeren die als invariante maat onze doeldistributie hebben. Hoe kunnen wij de keten irreducibel en aperiodiek krijgen? Zijn deze eisen in het algemeen te sterk? Verder zagen we in voorbeeld 3.8 dat de keuze van verschillende keuzes van σk tot een sneller resultaat konden leiden. Bestaat er een optimale σk ? Stel nou dat men teveel in een toestand blijft stilstaan in de Markov Chain. Dit wordt ook wel aangeduid met slow mixing. Is er een mogelijkheid dit soort momenten tegen te houden? Laat Fn (t) de benadering van de posterior verdeling zijn. Dan gegeven een ε > 0, kan ik dan een N ∈ N vinden zo dat voor alle n > N kFn (t) − Π(θ ≤ t|X n )k < ε? Een andere vraag is hoe snel de Gibbs Sampler convergeert vergeleken met het symmetrische Metropolis Hastings algoritme. Zijn er combinaties van beide algoritmes mogelijk, zo dat we een sneller een benadering van ons probleem kunnen vinden? Hoe wordt dat verschil tussen de benadering Fn en Π gemeten? Gebruikt men hiervoor de sup-norm of kan men een voordeligere norm definiëren? Bestaat er een bovengrens voor kFn − Π(θ ≤ t|X n )k, zoals dat bijvoorbeeld te vinden is in de numerieke wiskunde waar men de fout tussen een interpolatiepolynoom en een continue functie kan aangeven. Voor de Gibbs Sampler maken wij gebruik van conditionele verdelingen. Gegeven een posterior is het dan mogelijk deze zo in conditionele kansdichtheden te splitsen zo dat wij altijd de Gibbs Sampler kunnen gebruiken? Samenvattend kan men zeggen dat deze scriptie pas het topje van de ijsberg is en men nog veel dieper in de materie moet gaan om dergelijke vragen nauwkeurig te kunnen beantwoorden.

40

Populaire Samenvatting In deze scriptie bespreken wij Monte Carlo Markov Chains met hun toepassingen in de Bayesiaanse statistiek. Monte Carlo Markov Chains worden gebruikt om met de computer kansverdelingen te simuleren. De eerste Monte Carlo methoden werden in de 40er jaren bedacht door John van Neumann, Stanislaw Ulam en Nicholas Metropolis. De methodes zijn genoemd naar het Monte Carlo Casino, waar de oom van Ulam met grote regelmaat zijn geld verloor [6].

Figuur 4.1: Een foto van Nicholas Metropolis. Hij is een van de bedenkers van Monte Carlo methoden. Na hem is onder andere het symmetrische Metropolis Hastings algoritme genoemd [7]. De Monte Carlo methoden hebben verschillende toepassingen. Als we voor bepaalde stochasten bijvoorbeeld de variantie of de verwaching niet kunnen bepalen, kunnen Monte Carlo methoden zoals Monte Carlo integratie uitkomst bieden. We kunnen echter nog veel meer dan dat. In deze scriptie hebben we enkele simulatiemethoden zoals inversie en het verwerpingsalgoritme behandeld. Deze waren echter niet voldoende om complexe kansexperimenten zoals het benaderen van posteriors te simuleren. We hebben daarom de theorie van Markov ketens bestudeerd. Markov ketens beschrijven stochastische processen, die tussen verschillende toestanden springen. Deze sprongen 41

zijn geheugenloos. Dat wil zeggen dat een sprong naar een andere toestand niet afhangt van in welke toestand(en) men eerder verkeerde. De reden om naar Markov ketens te kijken, is dat Markov ketens een invariante maat hebben waar de keten als hij irreducibel en aperiodiek is naar toe convergeert. We construeren dan Markov ketens die als invariante maat de posterior verdeling hebben. Deze ketens noemen wij Monte Carlo Markov Chains. Tot slot geven wij een algoritme van een Monte Carlo Markov Chain om ons inleidende probleem 1.1 op te lossen. We kiezen als MCMC het symmetrische Metropolis Hastings algoritme: f u n c t i o n [ v a r ] = s y m m e t r o p o l i s h a s t i n g s m e a n s t u d e n t ( nu , N, sigma ) t h e t a ( 1 ) =0; f o r i =2:N Y=normrnd ( t h e t a ( i −1) , sigma ) ; % p r o p o s a l v o o r Y a l p h a=min ( ((1+(Yˆ2/ nu ) ) /(1+( t h e t a ( i −1) ˆ2/ nu ) ) ) ˆ( −(( nu+1) / 2 ) ) , 1 ) ; u=rand ; i f u
Figuur 4.2: Het symmetrische Metropolis Hastings voor de Student PN algoritme 1 2 t-verdeling in MATLAB. We schatten met N i=1 θi de variantie. De variantie van de Student t-verdeling is voor ν > 2 bekend en wordt ν gegeven door ν−2 . We kiezen voor onze MATLAB functie verschillende ν en herhalen het algoritme N = 10000 keer. Met de bekendheid van de variantie kunnen we inschatten hoe goed onze benadering is. We geven dit weer in de volgende tabel: ν

N

schatting variantie

variantie

3

10000

3.2425

3

4

10000

2.0589

2

5

10000

1.6576

1.666...

6

10000

1.6018

1.5

10

10000

1.2712

1.25

20

10000

1.0836

1.111...

42

Bibliografie [1] Christian P. Robert, The Bayesian Choice, 2nd Edition, Springer, Paris (2007) [2] Gareth O. Roberts, Jeffrey Rosenthal, General state space Markov chains and MCMC algorithms, Probability Surveys, Vol.1 [3] Bert van Es, Monte Carlo Markov Chain Simulation, 2011 [4] Bas Kleijn, Syllabus Bayesian Statistics, collegedictaat, UvA, Amsterdam (2011) [5] James R. Norris, Markov Chains, Cambridge Series in Statistical and Probabilistic Mathematics, Cambridge University Press, Cambridge (1998) [6] http://en.wikipedia.org/wiki/Monte_Carlo_method [7] http://www.lanl.gov/asc/metropolis.shtml

43

Monte Carlo Markov Chains voor Bayesiaanse statistiek

Recommend Documents