Thuis
Contacten

    Hoofdpagina


Docent: drs. Rob Flohr

Dovnload 2 Mb.

Docent: drs. Rob Flohr



Pagina8/25
Datum05.12.2018
Grootte2 Mb.

Dovnload 2 Mb.
1   ...   4   5   6   7   8   9   10   11   ...   25

Markov Chain Monte Carlo (MCMC) simulaties:
Omdat directe trekking uit de 'target distribution' (de posterior kansverdeling) soms niet mogelijk is worden de trekkingen gekoppeld: elke trekking is (alleen) afhankelijk van de voorgaande trekking, dit levert een Markov keten op (vgl. de ' random walk' zoals de wandeling van de dronkaard).
Voordeel: onder bepaalde voorwaarden convergeert een Markov keten naar een stabiele verdeling die de gewenste posterior verdeling benadert.
Er bestaan verschillende algoritmen om zo'n MCMC-simulatie uit te voeren, wij zullen de binnen de Bayesiaanse statistiek dominante Gibbs sampler gebruiken (is een variant op het zogeheten Metropolis-Hastings aloritme). De Gibbs sampler (vandaar de benaming BUGS: Bayesian Inference Using the Gibbs Sampler) kan gebruikt worden wanneer de posterior verdeling volledig bepaald wordt door de betreffende conditionele kansverdelingen.
We zullen tijdens deze cursus niet verder op de theoretische aspecten van de Gibbs sampler ingaan.

Er zijn verschillende softwaremogelijkheden om MCMC-simulaties uit te voeren, wij zullen ons beperken tot R-packages en het programma WinBUGS (zijn momenteel beide op de computers van de NHL geïnstalleerd).


VOORBEELDEN VAN TOEPASSINGEN VAN MCMC-pack:

A)
Normaal verdeelde prior en normaal verdeelde data (het voorbeeld onder punt 2) hierboven, p. 5-7)

> local({pkg <- select.list(sort(.packages(all.available = TRUE)),graphics=TRUE)

+ if(nchar(pkg)) library(pkg, character.only=TRUE)})

Loading required package: coda

Loading required package: lattice

Loading required package: MASS

##

## Markov Chain Monte Carlo Package (MCMCpack)



## Copyright (C) 2003-2014 Andrew D. Martin, Kevin M. Quinn, and Jong Hee Park

##

## Support provided by the U.S. National Science Foundation



## (Grants SES-0350646 and SES-0350613)

##

Warning messages:



1: package ‘MCMCpack’ was built under R version 2.15.3

2: package ‘coda’ was built under R version 2.15.3

3: package ‘lattice’ was built under R version 2.15.3

> library(MCMCpack)

> help(MCnormalnormal)

starting httpd help server ... done



de data: 100 aselecte trekkingen uit een normale verdeling met gemiddelde 68 en stand.dev. 4

> y=rnorm(100,68,4)

> y

[1] 72.45049 70.13388 76.40175 60.00040 63.67043 68.63739 65.72268 63.19363



[9] 63.38458 71.01923 65.75988 61.66443 69.36793 69.23036 61.45558 59.52161

[17] 73.48183 69.75289 73.14387 69.62487 64.32003 75.41437 68.04831 68.89889

[25] 62.40103 63.48986 62.90070 64.15750 73.80653 65.27796 69.30066 70.63795

[33] 66.22572 61.78972 75.94691 66.33717 71.93241 73.59413 69.99495 64.73935

[41] 68.98661 71.62322 65.54186 74.44156 62.25272 67.09229 64.00076 65.88842

[49] 70.26691 70.82734 68.76733 71.07976 63.92203 66.85162 62.97474 67.22301

[57] 70.25451 70.23980 59.98572 63.63328 73.75508 65.36202 64.97579 72.85769

[65] 73.65736 68.91184 66.09737 70.67694 71.79478 67.31302 71.03973 65.89309

[73] 68.72657 62.86158 64.70197 58.94994 74.83673 63.20605 61.81223 73.14511

[81] 69.13543 69.24244 64.64135 71.32216 67.98028 73.20444 65.53164 65.74803

[89] 68.63386 59.28228 69.85935 69.36069 66.71530 62.37406 64.12677 67.76551

[97] 65.47879 72.60752 70.69415 69.34366

> MCnormalnormal(y,16,72,36,5000)

y = de data
16 = de (bekende) variantie van de data y, dus 4^2=16
72 = de prior mean (het gemiddelde van de normaal verdeelde prior)
36 = de variantie van de normaal verdeelde prior
5000 = het aantal MCMC-trekkingen

Markov Chain Monte Carlo (MCMC) output:

Start = 1

End = 5000

Thinning interval = 1

mu

[1,] 68.00321



[2,] 67.59034

[3,] 68.04725

[4,] 66.94642

[5,] 67.97339

[6,] 68.04981

[7,] 67.87571

[8,] 67.62356

[9,] 67.32878

[10,] 67.25697

[11,] 68.05200

[12,] 67.04024

[13,] 68.01235

[14,] 67.23414

[15,] 67.61986

[16,] 67.84259

[17,] 67.16493

[18,] 67.86425

[19,] 67.42491

[20,] 67.27779
.
.
.
.

[4994,] 67.80474

[4995,] 67.83989

[4996,] 67.04659

[4997,] 68.28962

[4998,] 67.41418

[4999,] 67.09619

[5000,] 68.27077

attr(,"title")

[1] "MCnormalnormal Posterior Sample"

> summary(MCnormalnormal(y,16,72,36,5000))
Iterations = 1:5000

Thinning interval = 1

Number of chains = 1

Sample size per chain = 5000


1. Empirical mean and standard deviation for each variable,

plus standard error of the mean:


Mean SD Naive SE Time-series SE

67.673629 0.398766 0.005639 0.005639


2. Quantiles for each variable:
2.5% 25% 50% 75% 97.5%

66.89 67.40 67.68 67.94 68.45

B)
Voorbeeld Huiswerkopgave over Harry Potter m.b.v. MCMC-simulaties:

Met MCMCpack:


> local({pkg <- select.list(sort(.packages(all.available = TRUE)),graphics=TRUE)

+ if(nchar(pkg)) library(pkg, character.only=TRUE)})

Loading required package: coda

Loading required package: lattice

Loading required package: MASS

##

## Markov Chain Monte Carlo Package (MCMCpack)



## Copyright (C) 2003-2014 Andrew D. Martin, Kevin M. Quinn, and Jong Hee Park

##

## Support provided by the U.S. National Science Foundation



## (Grants SES-0350646 and SES-0350613)

##

Warning messages:



1: package ‘MCMCpack’ was built under R version 2.15.3

2: package ‘coda’ was built under R version 2.15.3

3: package ‘lattice’ was built under R version 2.15.3

> library(MCMCpack)

> harrypotter=MCbinomialbeta(y=400,n=1000,alpha=1,beta=1,mc=10000)

> summary(harrypotter)


Iterations = 1:10000

Thinning interval = 1

Number of chains = 1

Sample size per chain = 10000


1. Empirical mean and standard deviation for each variable,

plus standard error of the mean:


Mean SD Naive SE Time-series SE

0.4001817 0.0156443 0.0001564 0.0001564


2. Quantiles for each variable:
2.5% 25% 50% 75% 97.5%

0.3694 0.3896 0.4001 0.4108 0.4312

ad 4) installatie plus instructie R en gebruik van R-packages
(o.a.:
-install package -> load package
-library (naam package)
-library(help=naam package)
-help(naam specifieke R-opdracht)

(voorbeeld MCMCpack)

Huiswerkopgave voor les 3

I) Bestudeer hoofdstuk 4 uit De Bayesiaanse benadering

II) Bekijk het onderstaande vraagstuk.

Het wekelijkse aantal verkeersongevallen op de snelweg A.. tussen de plaatsen B. en C. volgt een Poisson-verdeling. Vier studenten gaan de komende acht weken het aantal verkeersongevallen op dat gedeelte van die snelweg vaststellen en de data op een Bayesiaanse manier analyseren.


Op basis van gegevens uit het verleden gaan ze uit van een gemiddeld aantal verkeersongevallen van 2.5 per week met een standaarddeviatie van 1. Ze drukken deze 'prior belief' uit middels een gamma-verdeling.
Nadat de acht weken zijn verstreken blijkt het aantal verkeersongevallen per week als volgt te zijn:

3, 2, 0, 8, 2, 4, 6, 1

Geef de karakteristieken van de betreffende posterior verdeling voor het aantal verkeersongevallen op deze snelweg
a)
langs analytische/theoretische weg, op basis van conjugacy, zonder software
b)
langs analytische/theoretische weg, op basis van conjugacy, maar wel met gebruik van relevante software
c)
met behulp van MCMC simulaties
d)
vergelijk de uitkomsten

Aanwijzing:


Om de parameters (, ) voor een gamma-prior te vinden kun je de volgende formules gebruiken:


Uitwerking:



a)

Aanwijzing:een prior met successen ( = aantal successen) geeft een posterior .
In het voorbeeld:
gamma(6.25, 2.5) gecombineerd met de data(successen) (3, 2, 0, 8, 2, 4, 6, 1)
geeft een gamma posterior ((6.25 + som data=6.25+26=32.25), ( 2.5 + aantal successen = 2.5+8=10.5)), dus gamma(32.25,10.5).
Gemiddelde: 32.25/10.5 = 3.07143, variantie = 32.25 / (10.5)^2 = 32.25/110.25 = 0.292517 -> stand. deviatie = 0.540848

Credible intervals met R:


99% credible interval voor mu: 1.9546 / 4.4676
>qgamma(0.01,32.25,10.5)
1.954566
>qgamma(0.99,32.25,10.5)
4.467603

95% credible interval voor mu: 2.1043 / 4.2186


>qgamma(0.025,32.25,10.5)
2.104253
>qgamma(0.975,32.25,10.5)
4.218608
b)

> local({pkg <- select.list(sort(.packages(all.available = TRUE)),graphics=TRUE)

+ if(nchar(pkg)) library(pkg, character.only=TRUE)})
Attaching package: ‘Bolstad’
The following object(s) are masked _by_ ‘.GlobalEnv’:
binodp
Warning message:

package ‘Bolstad’ was built under R version 2.15.3

> library(Bolstad)

> help(poisgamp)

starting httpd help server ... done

> y=c(3,2,0,8,2,4,6,1)

> results=poisgamp(y,6.25,2.5,plot=T)

Summary statistics for data

---------------------------

Number of observations: 8

Sum of observations: 26
Summary statistics for posterior

--------------------------------

Shape parameter r: 32.25

Rate parameter v: 10.5

99% credible interval for mu: [ 1.86 , 4.64 ]

> # mean=32.25/10.5 = 3.0714



> # variantie = 32.25/(10.5)^2 = 0.2925170068 -> sd = 0.5408

-----------------------------------------------------------------------------------------------------------------------------------

c)

> local({pkg <- select.list(sort(.packages(all.available = TRUE)),graphics=TRUE)



+ if(nchar(pkg)) library(pkg, character.only=TRUE)})

Loading required package: coda

Loading required package: lattice

Loading required package: MASS

##

## Markov Chain Monte Carlo Package (MCMCpack)



## Copyright (C) 2003-2014 Andrew D. Martin, Kevin M. Quinn, and Jong Hee Park

##

## Support provided by the U.S. National Science Foundation



## (Grants SES-0350646 and SES-0350613)

##

Warning messages:



1: package ‘MCMCpack’ was built under R version 2.15.3

2: package ‘coda’ was built under R version 2.15.3

3: package ‘lattice’ was built under R version 2.15.3

> library(MCMCpack)

> verkeersongevallen=c(3,2,0,8,2,4,6,1)

> posterior=MCpoissongamma(verkeersongevallen,6.25,2.5,5000)

> summary(posterior)
Iterations = 1:5000

Thinning interval = 1

Number of chains = 1

Sample size per chain = 5000


1. Empirical mean and standard deviation for each variable,

plus standard error of the mean:


Mean SD Naive SE Time-series SE

3.065579 0.535144 0.007568 0.007568


2. Quantiles for each variable:
2.5% 25% 50% 75% 97.5%

2.092 2.691 3.031 3.410 4.185


> plot(posterior)

>



Samengevat:


Op basis van gegevens uit het verleden (gemiddeld aantal verkeersongevallen van 2.5 per week met een standaarddeviatie van 1 ) en van recente data ( 3, 2, ....., 6, 1) kunnen we afleiden dat de kans 95% is dat het aantal wekelijkse ongevallen zal liggen tussen 2.1 en 4.2.

In het onderstaande gaan we uit van een lambda van 3.07.


a) In de week na het onderzoek doet zich slechts één ongeval voor. Is dit uitzonderlijk te noemen?
b) In de periode daarna is het aantal ongevallen 4 of meer. is dit uitzonderlijk te noemen?

Antwoorden:


a)
>dpois(1,3.07)
0.1425129

Nee, het kan aan toeval worden toegeschreven.

b)
>1-ppois(4,3.07)
0.196633

Nee, het kan aan toeval worden toegeschreven.


Les 3


Onderwerpen:

1)
a) groepsoefening: uitleg significantiebegrip aan een leek op statistisch gebied

b) belang van Bayesiaanse statistiek

c) andere R-packages voor Bayesiaanse statistiek, bijv. LearnBayes, bayess, etc. (zie www.r-project.org, onder CRAN)


2) Statistiek in context (voordat we volgende week voluit met software - WinBUGS - aan de gang gaan)


1   ...   4   5   6   7   8   9   10   11   ...   25


Dovnload 2 Mb.