PCM — Oct 15, 2013, 9:11 PM
#-----------------------------------------------------------------------------
# Should I believe in Santa Claus? (Claus Moebus, 2013-10-15)
library(BRugs)
Welcome to BRugs connected to OpenBUGS version 3.2.2
#-----------------------------------------------------------------------------
# THE MODEL.
# Specify the model in BUGS language, but save it as a string in R:
modelString = "
# BUGS model specification begins ...
model {
# Prior Belief Distribution: No idea about the probability of his existence
theta ~ dbeta( priorHeads , priorTails )
priorHeads <- 1
priorTails <- 1
# Likelihood: If he exists with probability theta, then anyone
# I ask has the chance to see him with that probability theta
for ( i in 1:N ) {
seen[i] ~ dbern( theta )
}
# Posteriori Belief in Favor (flipping 'mental belief coin') P(0.5 <= theta)
mbCoin <- step(theta - 0.5) # E(mbCoin) = P(0.5 <= theta)
}
# ... BUGS model specification ends.
" # close quote to end modelString
# Write the modelString to a file, using R commands:
writeLines(modelString,con="model.txt")
# Use BRugs to send the model.txt file to BUGS, which checks the model syntax:
modelCheck( "model.txt" )
model is syntactically correct
#------------------------------------------------------------------------------
# THE DATA.
# Specify the data in R, using a list format compatible with BUGS:
# Each person in the sample with size N answers the question
# "Have you ever seen the true Santa Claus?" truthfully
dataList = list(
N = 5 ,
seen = c( 0, 0, 1, 0, 0 )
# N = 10 ,
# seen = c( 0, 0, 1, 0, 0, 0, 1, 0, 0, 0 )
)
# Use BRugs commands to put the data into a file and ship the file to BUGS:
modelData( bugsData( dataList ) )
data loaded
#------------------------------------------------------------------------------
# INTIALIZE THE CHAIN.
modelCompile() # BRugs command tells BUGS to compile the model.
model compiled
modelGenInits() # BRugs command tells BUGS to randomly initialize a chain.
initial values generated, model initialized
#------------------------------------------------------------------------------
# RUN THE CHAINS.
modelUpdate(1000) # burn in
1000 updates took 0 s
# BRugs tells BUGS to keep a record of the sampled "theta" values:
samplesSet( c("theta", "mbCoin"))
monitor set for variable 'theta'
monitor set for variable 'mbCoin'
# R command defines a new variable that specifies an arbitrary chain length:
chainLength = 10000
# BRugs tells BUGS to generate a MCMC chain:
modelUpdate( chainLength )
10000 updates took 0 s
#------------------------------------------------------------------------------
# EXAMINE THE RESULTS.
thetaSample = samplesSample( "theta") # BRugs asks BUGS for the sample values.
mbCoinSample = samplesSample( "mbCoin") # BRugs asks BUGS for the sample values.
thetaSummary = samplesStats( c("theta", "mbCoin")) # BRugs asks BUGS for summary statistics.
samplesStats("*")
mean sd MC_error val2.5pc median val97.5pc start sample
mbCoin 0.1087 0.3113 0.002827 0.0000 0.0000 1.0000 1001 10000
theta 0.2855 0.1607 0.001518 0.0431 0.2658 0.6452 1001 10000
# Make a graph using R commands:
windows(10,6)
layout( matrix( c(1,2) , nrow=1 ) )
plot( thetaSample[1:500] , 1:length(thetaSample[1:500]) , type="o" ,
xlim=c(0,1) , xlab=bquote(theta) , ylab="Position in Chain" ,
cex.lab=1.25 , main="BUGS Results", col="deepskyblue" )
plot( mbCoinSample[1:500] , 1:length(mbCoinSample[1:500]) , type="o" ,
xlim=c(0,1) , xlab=bquote(mbCoin) , ylab="Position in Chain" ,
cex.lab=1.25 , main="BUGS Results", col="deepskyblue" )
windows(10,6)
layout( matrix( c(1,2) , nrow=1 ) )
plotDensity("theta")
plotDensity("mbCoin")
source("plotPost.R")
plotPost( thetaSample , xlim=c(0,1), xlab="theta", main="theta-posterior" )
mean median mode hdiMass hdiLow hdiHigh compVal pcGTcompVal
theta 0.2855 0.2657 0.1751 0.95 0.01651 0.5908 NA NA
ROPElow ROPEhigh pcInROPE
theta NA NA NA
#-----------------------------------------------------------------------------
