Thuis
Contacten

    Hoofdpagina


Docent: drs. Rob Flohr

Dovnload 2 Mb.

Docent: drs. Rob Flohr



Pagina14/25
Datum05.12.2018
Grootte2 Mb.

Dovnload 2 Mb.
1   ...   10   11   12   13   14   15   16   17   ...   25

MODEL {

p ~ dbeta(a,b)

x ~ dbin(p,n)

}

DATA list(a=1,b=1,x=3140,n=13819)

Prior voor p


Likelihood/measurement model voor x




model is syntactically correct

data loaded

model compiled

initial values generated, model initialized
Dynamic trace
Kernel density
Node statistics

node mean sd MC error 2.5% median 97.5% start sample

p 0.2273 0.003596 1.184E-4 0.2206 0.2272 0.2343 1 1000






Ter vergelijking:

KLASSIEKE (FREQUENTISTISCHE) OPLOSSING:

De steekproefproportie is



Om het 95% betrouwbaarheidsinterval te vinden bepalen we eerst de standaardfout:



en vinden dan



hetgeen betekent dat we op basis van deze steekproef met een betrouwbaarheid van 95% schatten dat tussen de 22.0% en 23.4% van de betreffende studentenpopulatie 'frequent binge drinkers' zijn.

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

M.B.V. R:


> binom.test(3140,13819,0.22)

OUTPUT:


Exact binomial test
data: 3140 and 13819

number of successes = 3140, number of trials = 13819, p-value = 0.04102

alternative hypothesis: true probability of success is not equal to 0.22

95 percent confidence interval:

0.2202586 0.2343032

sample estimates:

probability of success

0.2272234


_______________________________________________________________


BAYESIAANSE OPLOSSING M.B.V. 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)

> posterior=MCbinomialbeta(3140,13819,alpha=1,beta=1,mc=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

2.272e-01 3.574e-03 5.054e-05 5.054e-05


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

0.2202 0.2248 0.2272 0.2296 0.2342


> plot(posterior)

>

(2)



Sla het odc-bestand 'cholesterol level children aged 2-14' op, open het in WinBUGS en laat het model draaien.
Aanwijzingen:
- laat de software de startwaarden genereren (gen inits)
- kies bij inference/samples - node de volgende parameters: population.mean, tau en sd
Output:

model is syntactically correct

data loaded

model compiled

initial values generated, model initialized

Time series


Kernel density
Node statistics

node mean sd MC error 2.5% median 97.5% start sample

pop.mean 207.1 3.201 0.04408 201.4 207.1 213.1 1 5000

sd 29.9 2.322 0.03498 26.04 29.77 34.47 1 5000

tau 0.001137 1.629E-4 2.4E-6 8.417E-4 0.001128 0.001475 1 5000

Les 5

De vier componenten van een Bayesiaanse analyse:

* een prior verdeling voor de parameter(s) (een parameter is een onbekende populatiegrootheid die geschat moet worden)
* de data
* een model dat de parameter(s) in verband brengt met/ relateert aan de data
* een posterior verdeling voor de parameter(s)

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


Voorbeelden met WinBUGS:

Een bioloog is geïnteresseerd in het gewicht (gemiddelde, spreiding) van een populatie slechtvalken. Er wordt een steekproef genomen van 10 slechtvalken, van elk exemplaar wordt het gewicht in grammen gemeten.


Welke uitspraak over het gewicht van de populatie kan op grond daarvan worden geformuleerd?

1. Klassiek: .........................


2. Bayesiaans:

# Inferring mean and standard deviation of a population distribution of peregrine falcons from a sample of observed independent data, assuming data following a Gaussian distribution.
list(x=c(657.33,593.92,634.68,589.31,647.24,582.21,571.38,623.63,561.90,586.39),n=10)
model{

for (i in 1:n) {

x[i]~dnorm(mu,lambda)

}
# Priors


mu~dunif(0,5000)

sigma~dunif(0,100)

lambda<- 1/pow(sigma,2)

}
WinBUGS OUTPUT:

model is syntactically correct

data loaded

model compiled

initial values generated, model initialized


Time series



Kernel density







Node statistics



node mean sd MC error 2.5% median 97.5% start sample

mu 604.7 12.75 0.1991 579.5 604.6 629.9 501 4500

sigma 38.97 11.06 0.239 23.59 37.04 66.7 501 4500
Klassiek:
Theorie/Formule:
unknown: one-sample test (one-sample statistic)



R-command:


t.test(data, mu= .... , conf.level=....)

> gewicht=c(657.33,593.92,634.68,589.31,647.24,582.21,571.38,623.63,561.90,586.39)

> t.test(gewicht,mu=400)


One Sample t-test
data: gewicht

t = 19.4554, df = 9, p-value = 1.158e-08

alternative hypothesis: true mean is not equal to 400

95 percent confidence interval:

580.9862 628.6118

sample estimates:

mean of x

604.799

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

vergelijken van twee behandelingen voor een bepaalde ziekte (ziekte is teruggekomen:return rates)

#Comparing results of two different treatments for a certain disease
list(k1=2,n1=23,k2=3,n2=18)
model{

k1~dbin(theta1,n1)

k2~dbin(theta2,n2)
#priors
theta1~dbeta(1,1)

theta2~dbeta(1,1)


#Difference return rates

delta<-theta1-theta2

}


WinBUGS OUTPUT:

model is syntactically correct

data loaded

model compiled

initial values generated, model initialized



Time series

Kernel density




Node statistics



node mean sd MC error 2.5% median 97.5% start sample

delta -0.07935 0.1094 0.001488 -0.3067 -0.07578 0.127 1 5000


Klassiek:

> prop.test(x=c(2,3),n=c(23,18),correct=FALSE)
2-sample test for equality of proportions without continuity

correction


data: c(2, 3) out of c(23, 18)

X-squared = 0.5992, df = 1, p-value = 0.4389

alternative hypothesis: two.sided

95 percent confidence interval:

-0.2868368 0.1274165

sample estimates:

prop 1 prop 2

0.08695652 0.16666667


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


THEORIE:

Statistiek: gericht op het kwantificeren van onzekerheid. Daarvoor hebben we een kansverdeling nodig.

Klassieke manier:
Het concept van de 'likelihood' werd in 1922 geïntroduceerd door Fisher. Dit drukt uit hoe plausibel de data zijn als functie van de parameter(s) van een statistisch (stochastisch) model. Als functie van de parameter krijgen we de 'likelihood function'.

Voorbeeld: 20 keer gooien van een munt, resultaat is 5 keer kop, dit zijn de data.

Parameter () is in dit voorbeeld : de proportie kop van de munt (d.w.z. de echte proportie kop van de munt, dus zoals die in werkelijkheid is en die we helaas nooit exact zullen weten)

Het statistisch model (dat verondersteld wordt de data te genereren) is in dit voorbeeld het binomiale model.

Als we op de horizontale as de mogelijke waarden van de parameter zetten en vervolgens op de verticale as de (op basis van het binomiale model berekende kansen op de data, d.w.z. 5 keer kop uit 20 worpen), dan krijg je de grafiek van de binomiale likelihood functie.
Het zal geen verbazing wekken dat de modus ligt bij .

De waarde van die de likelihood functie maximaliseert wordt de maximum likelihood estimate (MLE) genoemd en wordt doorgaans aangeduid door .

Het principe is:
Alle informatie ('all evidence') die we hebben verkregen uit ons onderzoek/experiment omtrent een onbekende parameter ligt besloten in de likelihood functie van , gegeven de data.

Vervolgens kunnen we een 95% betrouwbaarheidsinterval berekenen voor op basis van een steekproevenverdeling ('sampling distribution'). Deze steekproevenverdeling kan langs theoretische weg (Centrale Limiet Stelling) of langs empirische weg ('Bootstrap') gevonden worden. Uitgangspunt voor de schatting/betrouwbaarheidsinterval is de MLE, hier 0.25. De steekproevenverdeling is een kansverdeling op basis waarvan we onzekerheid kwantificeren.

Een voorspelling omtrent toekomstige data gebeurt op basis van de MLE en een 95% - 'predictive interval' dat 95% van de toekomstige data zal bevatten

Bayesiaanse manier:


Alle informatie waarover we beschikken, gebaseerd op prior informatie en op de data, ligt besloten in de posterior kansverdeling. Dit is de kansverdeling op basis waarvan we onzekerheid kwantificeren.

Op basis van de gevonden posterior distribution zetten we nu de stap naar de




Posterior predictive distribution

Posterior verdeling: deze bevat alle informatie over de parameter op basis van de prior en de data; we willen al deze informatie gebruiken om toekomstige data (afkomstig uit dezelfde populatie als waaruit de geobserveerde data afkomstig zijn) te voorspellen.

We beginnen met een discrete posterior verdeling, bijv. (zie vb munt) met 3 parameterwaarden met bijbehorende kansen, theta1, theta2, theta3 (theta=proportie kop)

We nemen theta1 en bepalen een kansverdeling voor het aantal successen, bijv. het gooien van 'kop', (gegeven een bepaald aantal worpen en het betreffende statistisch model, hier het binomiale model) op basis van theta1, vervolgens doen we hetzelfde voor theta2 en theta3 . Zo krijgen we drie verschillende kansverdelingen voor het aantal successen die we vervolgens sommeren, met in achtneming van de kansen op theta1, theta2 en theta3 (gewogen som) om de verdeling van te schatten:

Voor een continue posterior krijgen we dan:








In het voorbeeld van de munt vonden we (op basis van een uniforme prior) de volgende posterior verdeling:


De Bayesiaanse manier om te voorspellen gaat als volgt:

Wat we in feite doen is een groot aantal parameterwaarden aselect uit de posterior verdeling trekken en dan elke waarde gebruiken om het aantal successen (k) te bepalen op basis van ons binomiaal model (gerepliceerde waarden, 'replication').

De posterior predictive distribution kan gezien worden als een verzameling van voorspellingen ten aanzien van welke data op basis van het model verwacht kunnen worden.

Het concept van de posterior predictive distribution is ook relevant in de context van de vraag hoe adequaat een statistisch model is in het beschrijven van de werkelijkheid.

Wiskundig gezien wordt de posterior predictive distribution in dit voorbeeld gevonden door:

(voor alle waarden van )

Merk op dat we ons in deze benadering niet beperken tot de MLE () maar alle mogelijke waarden van de parameter bij onze voorspelling gebruiken.


Voorbeeld van de munt:

De posterior predictive distribution voor n=20 (gerepliceerde waarden voor k) levert het volgende op:

list(k=5,n=20)

model{

k~dbin(theta,n)



theta~dbeta(1,1)

postpredk~dbin(theta,n)


}

model is syntactically correct

data loaded

model compiled



initial values generated, model initialized
Time series



Kernel density



Node statistics


1   ...   10   11   12   13   14   15   16   17   ...   25

  • node mean sd MC error 2.5% median 97.5% start sample
  • > binom.test(3140,13819,0.22) OUTPUT
  • > library(MCMCpack)

  • Dovnload 2 Mb.