forked from boboppie/kruschke-doing_bayesian_data_analysis
-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy pathBernTwoFurrowsBugs.R
86 lines (75 loc) · 2.87 KB
/
BernTwoFurrowsBugs.R
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
graphics.off()
rm(list=ls(all=TRUE))
library(BRugs) # Kruschke, J. K. (2010). Doing Bayesian data analysis:
# A Tutorial with R and BUGS. Academic Press / Elsevier.
#------------------------------------------------------------------------------
# THE MODEL.
modelstring = "
# BUGS model specification begins here...
model {
# Likelihood. Each flip is Bernoulli.
for ( i in 1 : N1 ) { y1[i] ~ dbern( theta1 ) }
for ( i in 1 : N2 ) { y2[i] ~ dbern( theta2 ) }
# Prior. Curved scallops!
x ~ dunif(0,1)
y ~ dunif(0,1)
N <- 4
xt <- sin( 2*3.141593*N * x ) / (2*3.141593*N) + x
yt <- 3 * y + (1/3)
xtt <- pow( xt , yt )
theta1 <- xtt
theta2 <- y
}
# ... end BUGS model specification
" # close quote for modelstring
# Write model to a file:
.temp = file("model.txt","w") ; writeLines(modelstring,con=.temp) ; close(.temp)
# Load model file into BRugs and check its syntax:
modelCheck( "model.txt" )
#------------------------------------------------------------------------------
# THE DATA.
# Specify the data in a form that is compatible with BRugs model, as a list:
datalist = list(
N1 = 7 ,
y1 = c( 1,1,1,1,1,0,0 ) ,
N2 = 7 ,
y2 = c( 1,1,0,0,0,0,0 )
)
# Get the data into BRugs:
modelData( bugsData( datalist ) )
#------------------------------------------------------------------------------
# INTIALIZE THE CHAIN.
modelCompile()
modelGenInits()
#------------------------------------------------------------------------------
# RUN THE CHAINS.
samplesSet( c( "theta1" , "theta2" ) ) # Keep a record of sampled "theta" values
chainlength = 10000 # Arbitrary length of chain to generate.
modelUpdate( chainlength ) # Actually generate the chain.
#------------------------------------------------------------------------------
# EXAMINE THE RESULTS.
theta1Sample = samplesSample( "theta1" ) # Put sampled values in a vector.
theta2Sample = samplesSample( "theta2" ) # Put sampled values in a vector.
source("plotChains.R")
plotChains("theta1")
plotChains("theta2")
# Plot the trajectory of the last sampled values.
windows()
par( pty="s" )
nToPlot=2000
plot( theta1Sample[(chainlength-nToPlot+1):chainlength] ,
theta2Sample[(chainlength-nToPlot+1):chainlength] , type = "p" ,
xlim = c(0,1) , xlab = bquote(theta[1]) , ylim = c(0,1) ,
ylab = bquote(theta[2]) , main="BUGS Result" )
# Display means in plot.
theta1mean = mean(theta1Sample)
theta2mean = mean(theta2Sample)
if (theta1mean > .5) { xpos = 0.0 ; xadj = 0.0
} else { xpos = 1.0 ; xadj = 1.0 }
if (theta2mean > .5) { ypos = 0.0 ; yadj = 0.0
} else { ypos = 1.0 ; yadj = 1.0 }
text( xpos , ypos ,
bquote(
"M=" * .(signif(theta1mean,3)) * "," * .(signif(theta2mean,3))
) ,adj=c(xadj,yadj) ,cex=1.5 )
dev.copy2eps(file="BernTwoFurrowsBugs.eps")