Skip to content

Commit

Permalink
updated codes
Browse files Browse the repository at this point in the history
  • Loading branch information
@boboppie
boboppie committed Feb 8, 2013
1 parent 4879c82 commit fa4c5ba
Show file tree
Hide file tree
Showing 15 changed files with 234 additions and 156 deletions.
25 changes: 12 additions & 13 deletions BayesUpdate.R
View file
Original file line number Diff line number Diff line change
Expand Up @@ -11,10 +11,10 @@ pTheta = pmin( Theta , 1-Theta ) # Makes a triangular belief distribution.
pTheta = pTheta / sum( pTheta ) # Makes sure that beliefs sum to 1.

# Specify the data. To produce the examples in the book, use either
# Data = c(1,1,1,0,0,0,0,0,0,0,0,0) or Data = c(1,0,0,0,0,0,0,0,0,0,0,0).
Data = c(1,1,1,0,0,0,0,0,0,0,0,0)
nHeads = sum( Data == 1 )
nTails = sum( Data == 0 )
# Data = c(rep(1,3),rep(0,9)) or Data = c(rep(1,1),rep(0,11)).
Data = c(rep(1,3),rep(0,9))
nHeads = sum( Data )
nTails = length( Data ) - nHeads

# Compute the likelihood of the data for each value of theta:
pDataGivenTheta = Theta^nHeads * (1-Theta)^nTails
Expand All @@ -24,7 +24,8 @@ pData = sum( pDataGivenTheta * pTheta )
pThetaGivenData = pDataGivenTheta * pTheta / pData # This is Bayes' rule!

# Plot the results.
windows(7,10) # create window of specified width,height inches.
source("openGraphSaveGraph.R") # read in graph functions
openGraph(width=7,height=10,mag=0.7) # open a window for the graph
layout( matrix( c( 1,2,3 ) ,nrow=3 ,ncol=1 ,byrow=FALSE ) ) # 3x1 panels
par(mar=c(3,3,1,0)) # number of margin lines: bottom,left,top,right
par(mgp=c(2,1,0)) # which margin lines to use for labels
Expand All @@ -34,30 +35,28 @@ par(mai=c(0.5,0.5,0.3,0.1)) # margin size in inches: bottom,left,top,right
plot( Theta , pTheta , type="h" , lwd=3 , main="Prior" ,
xlim=c(0,1) , xlab=bquote(theta) ,
ylim=c(0,1.1*max(pThetaGivenData)) , ylab=bquote(p(theta)) ,
cex.axis=1.2 , cex.lab=1.5 , cex.main=1.5 )
cex.axis=1.2 , cex.lab=1.5 , cex.main=1.5 , col="skyblue" )

# Plot the likelihood:
plot( Theta , pDataGivenTheta , type="h" , lwd=3 , main="Likelihood" ,
xlim=c(0,1) , xlab=bquote(theta) ,
ylim=c(0,1.1*max(pDataGivenTheta)) , ylab=bquote(paste("p(D|",theta,")")),
cex.axis=1.2 , cex.lab=1.5 , cex.main=1.5 )
cex.axis=1.2 , cex.lab=1.5 , cex.main=1.5 , col="skyblue" )
text( .55 , .85*max(pDataGivenTheta) , cex=2.0 ,
bquote( "D=" * .(nHeads) * "H," * .(nTails) * "T" ) , adj=c(0,.5) )

# Plot the posterior:
plot( Theta , pThetaGivenData , type="h" , lwd=3 , main="Posterior" ,
xlim=c(0,1) , xlab=bquote(theta) ,
ylim=c(0,1.1*max(pThetaGivenData)) , ylab=bquote(paste("p(",theta,"|D)")),
cex.axis=1.2 , cex.lab=1.5 , cex.main=1.5 )
cex.axis=1.2 , cex.lab=1.5 , cex.main=1.5 , col="skyblue" )
text( .55 , .85*max(pThetaGivenData) , cex=2.0 ,
bquote( "p(D)=" * .(signif(pData,3)) ) , adj=c(0,.5) )

# Save the plot as an EPS file.
# Save the plot
if ( nThetaVals == 3 ) { modeltype = "simpleModel" }
if ( nThetaVals == 63 ) { modeltype = "complexModel" }
if ( nHeads == 3 & nTails == 9 ) { datatype = "simpleData" }
if ( nHeads == 1 & nTails == 11 ) { datatype = "complexData" }
filename = paste( "BayesUpdate_" ,modeltype ,"_" ,datatype ,".eps" ,sep="" )
# The command dev.copy2eps, used below, doesn't work on all systems.
# Try help("dev.copy2eps") for info about saving graphs in other file formats.
dev.copy2eps( file=filename )
filename = paste( "BayesUpdate_" ,modeltype ,"_" ,datatype,sep="" )
saveGraph( file=filename , type="eps" )
20 changes: 11 additions & 9 deletions BernBeta.R
View file
Original file line number Diff line number Diff line change
@@ -1,4 +1,4 @@
BernBeta = function( priorShape , dataVec , credMass=0.95 , saveGraph=F ) {
BernBeta = function( priorShape , dataVec , credMass=0.95 , saveGr=FALSE ) {
# Bayesian updating for Bernoulli likelihood and beta prior.
# Input arguments:
# priorShape
Expand All @@ -18,6 +18,7 @@ BernBeta = function( priorShape , dataVec , credMass=0.95 , saveGraph=F ) {
# You will need to "source" this function before using it, so R knows
# that the function exists and how it is defined.


# Check for errors in input arguments:
if ( length(priorShape) != 2 ) {
stop("priorShape must have two components.") }
Expand Down Expand Up @@ -54,15 +55,16 @@ pDataGivenTheta = Theta^z * (1-Theta)^(N-z)
# Compute the posterior at each value of theta.
pThetaGivenData = dbeta( Theta , a+z , b+N-z )
# Open a window with three panels.
windows(7,10)
source("openGraphSaveGraph.R") # read in graph functions
openGraph(width=7,height=10,mag=0.7) # open a window for the graph
layout( matrix( c( 1,2,3 ) ,nrow=3 ,ncol=1 ,byrow=FALSE ) ) # 3x1 panels
par( mar=c(3,3,1,0) , mgp=c(2,1,0) , mai=c(0.5,0.5,0.3,0.1) ) # margin specs
maxY = max( c(pTheta,pThetaGivenData) ) # max y for plotting
# Plot the prior.
plot( Theta , pTheta , type="l" , lwd=3 ,
xlim=c(0,1) , ylim=c(0,maxY) , cex.axis=1.2 ,
xlab=bquote(theta) , ylab=bquote(p(theta)) , cex.lab=1.5 ,
main="Prior" , cex.main=1.5 )
main="Prior" , cex.main=1.5 , col="skyblue" )
if ( a > b ) { textx = 0 ; textadj = c(0,1) }
else { textx = 1 ; textadj = c(1,1) }
text( textx , 1.0*max(pThetaGivenData) ,
Expand All @@ -73,7 +75,7 @@ plot( Theta , pDataGivenTheta , type="l" , lwd=3 ,
xlim=c(0,1) , cex.axis=1.2 , xlab=bquote(theta) ,
ylim=c(0,1.1*max(pDataGivenTheta)) ,
ylab=bquote( "p(D|" * theta * ")" ) ,
cex.lab=1.5 , main="Likelihood" , cex.main=1.5 )
cex.lab=1.5 , main="Likelihood" , cex.main=1.5 , col="skyblue" )
if ( z > .5*N ) { textx = 0 ; textadj = c(0,1) }
else { textx = 1 ; textadj = c(1,1) }
text( textx , 1.0*max(pDataGivenTheta) , cex=2.0 ,
Expand All @@ -82,7 +84,7 @@ text( textx , 1.0*max(pDataGivenTheta) , cex=2.0 ,
plot( Theta , pThetaGivenData ,type="l" , lwd=3 ,
xlim=c(0,1) , ylim=c(0,maxY) , cex.axis=1.2 ,
xlab=bquote(theta) , ylab=bquote( "p(" * theta * "|D)" ) ,
cex.lab=1.5 , main="Posterior" , cex.main=1.5 )
cex.lab=1.5 , main="Posterior" , cex.main=1.5 , col="skyblue" )
if ( a+z > b+N-z ) { textx = 0 ; textadj = c(0,1) }
else { textx = 1 ; textadj = c(1,1) }
text( textx , 1.00*max(pThetaGivenData) , cex=2.0 ,
Expand All @@ -101,10 +103,10 @@ text( hpdLim[1] , hpdHt , bquote(.(round(hpdLim[1],3))) ,
adj=c(1.1,-0.1) , cex=1.2 )
text( hpdLim[2] , hpdHt , bquote(.(round(hpdLim[2],3))) ,
adj=c(-0.1,-0.1) , cex=1.2 )
# Construct file name for saved graph, and save the graph.
if ( saveGraph ) {
filename = paste( "BernBeta_",a,"_",b,"_",z,"_",N,".eps" ,sep="")
dev.copy2eps( file = filename )
# Construct filename for saved graph, and save the graph.
if ( saveGr ) {
filename = paste( "BernBeta_",a,"_",b,"_",z,"_",N ,sep="")
saveGraph( file = filename , type="eps" )
}
return( postShape )
} # end of function
20 changes: 9 additions & 11 deletions BernBetaJagsFull.R
View file
Original file line number Diff line number Diff line change
@@ -1,8 +1,6 @@
rm(list = ls())
graphics.off()
if ( .Platform$OS.type != "windows" ) {
windows <- function( ... ) X11( ... )
}
source("openGraphSaveGraph.R")

require(rjags) # Kruschke, J. K. (2011). Doing Bayesian Data Analysis:
# A Tutorial with R and BUGS. Academic Press / Elsevier.
Expand Down Expand Up @@ -79,14 +77,14 @@ mcmcChain = as.matrix( codaSamples )
thetaSample = mcmcChain

# Make a graph using R commands:
windows(10,6)
openGraph(width=10,height=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="JAGS Results" )
cex.lab=1.25 , main="JAGS Results" , col="skyblue" )
source("plotPost.R")
histInfo = plotPost( thetaSample , xlim=c(0,1) , col="skyblue" )
savePlot( file="BernBetaJagsFull.eps" , type="eps" )
histInfo = plotPost( thetaSample , xlim=c(0,1) , xlab=bquote(theta) )
saveGraph( file="BernBetaJagsFullPost" , type="eps" )

# Posterior prediction:
# For each step in the chain, use posterior theta to flip a coin:
Expand All @@ -99,16 +97,16 @@ for ( stepIdx in 1:chainLength ) {
# Jitter the 0,1 y values for plotting purposes:
yPredJittered = yPred + runif( length(yPred) , -.05 , +.05 )
# Now plot the jittered values:
windows(5,5.5)
openGraph(width=5,height=5.5)
plot( thetaSample[1:500] , yPredJittered[1:500] , xlim=c(0,1) ,
main="posterior predictive sample" ,
xlab=expression(theta) , ylab="y (jittered)" )
points( mean(thetaSample) , mean(yPred) , pch="+" , cex=2 )
xlab=expression(theta) , ylab="y (jittered)" , col="skyblue" )
points( mean(thetaSample) , mean(yPred) , pch="+" , cex=2 , col="skyblue" )
text( mean(thetaSample) , mean(yPred) ,
bquote( mean(y) == .(signif(mean(yPred),2)) ) ,
adj=c(1.2,.5) )
text( mean(thetaSample) , mean(yPred) , srt=90 ,
bquote( mean(theta) == .(signif(mean(thetaSample),2)) ) ,
adj=c(1.2,.5) )
abline( 0 , 1 , lty="dashed" , lwd=2 )
savePlot( file="BernBetaJagsPost.eps" , type="eps" )
saveGraph( file="BernBetaJagsFullPred" , type="eps" )
40 changes: 20 additions & 20 deletions BernBetaMuKappaJags.R
View file
Original file line number Diff line number Diff line change
@@ -1,9 +1,6 @@
rm(list = ls())
graphics.off()
if ( .Platform$OS.type != "windows" ) {
windows <- function( ... ) X11( ... )
}

source("openGraphSaveGraph.R")
require(rjags) # Kruschke, J. K. (2011). Doing Bayesian Data Analysis:
# A Tutorial with R and BUGS. Academic Press / Elsevier.
#------------------------------------------------------------------------------
Expand Down Expand Up @@ -106,11 +103,13 @@ codaSamples = coda.samples( jagsModel , variable.names=parameters ,
#------------------------------------------------------------------------------
# EXAMINE THE RESULTS.

checkConvergence = F
checkConvergence = FALSE
if ( checkConvergence ) {
show( summary( codaSamples ) )
plot( codaSamples , ask=T )
autocorr.plot( codaSamples , ask=T )
openGraph(width=7,height=7)
autocorr.plot( codaSamples[[1]] , ask=FALSE )
show( gelman.diag( codaSamples ) )
effectiveChainLength = effectiveSize( codaSamples )
show( effectiveChainLength )
}

# Convert coda-object codaSamples to matrix object for easier handling.
Expand All @@ -130,7 +129,7 @@ for ( coinIdx in 1:nCoins ) {
# Make a graph using R commands:
source("plotPost.R")
if ( nCoins <= 5 ) { # Only make this figure if there are not too many coins
windows(3.2*3,2.5*(1+nCoins))
openGraph(width=3.2*3,height=2.5*(1+nCoins))
layout( matrix( 1:(3*(nCoins+1)) , nrow=(nCoins+1) , byrow=T ) )
par(mar=c(2.95,2.95,1.0,0),mgp=c(1.35,0.35,0),oma=c( 0.1, 0.1, 0.1, 0.1))
nPtsToPlot = 500
Expand All @@ -139,12 +138,11 @@ if ( nCoins <= 5 ) { # Only make this figure if there are not too many coins
plot( muSample[plotIdx] , kappaSample[plotIdx] , type="p" , ylim=kPltLim ,
xlim=c(0,1) , xlab=expression(mu) , ylab=expression(kappa) , cex.lab=1.5 ,
col="skyblue" )
plotPost( muSample , xlab="mu" , xlim=c(0,1) , main="" , breaks=20 , col="skyblue")
plotPost( kappaSample , xlab="kappa" , main="" , breaks=20 , col="skyblue" ,
showMode=T )
plotPost( muSample , xlab=bquote(mu) , xlim=c(0,1) , main="")
plotPost( kappaSample , xlab=bquote(kappa) , main="" , showMode=TRUE )
for ( coinIdx in 1:nCoins ) {
plotPost( thetaSample[coinIdx,] , xlab=paste("theta",coinIdx,sep="") ,
xlim=c(0,1) , main="" , breaks=20 , col="skyblue")
plotPost( thetaSample[coinIdx,] , xlab=bquote(theta[.(coinIdx)]) ,
xlim=c(0,1) , main="" )
plot( thetaSample[coinIdx,plotIdx] , muSample[plotIdx] , type="p" ,
xlim=c(0,1) , ylim=c(0,1) , cex.lab=1.5 ,
xlab=bquote(theta[.(coinIdx)]) , ylab=expression(mu) , col="skyblue")
Expand All @@ -155,11 +153,13 @@ if ( nCoins <= 5 ) { # Only make this figure if there are not too many coins
} # end if ( nCoins <= ...

#
windows(width=7,height=5)
openGraph(width=7,height=5)
layout( matrix( 1:4 , nrow=2 , byrow=T ) )
par(mar=c(2.95,2.95,1.0,0),mgp=c(1.35,0.35,0),oma=c( 0.1, 0.1, 0.1, 0.1) )
plotPost( muSample , xlab="mu" , main="" , breaks=20 , compVal=0.5 , col="skyblue")
plotPost( kappaSample , xlab="kappa" , main="" , breaks=20 , HDItextPlace=.1 , col="skyblue")
plotPost( thetaSample[1,] , xlab="theta1" , main="" , breaks=20 , compVal=0.5 , col="skyblue")
# Uncomment the following if using therapeutic touch data:
#plotPost( thetaSample[28,] , xlab="theta28" , main="" , breaks=20 , compVal=0.5 , col="skyblue")
plotPost( muSample , xlab=bquote(mu) , main="" , compVal=0.5 )
plotPost( kappaSample , xlab=bquote(kappa) , main="" , HDItextPlace=.1 ,
showMode=TRUE)
plotPost( thetaSample[1,] , xlab=bquote(theta[1]) , main="" , compVal=0.5 )
lastIdx = length(z)
plotPost( thetaSample[lastIdx,] , xlab=bquote(theta[.(lastIdx)]) ,
main="" , compVal=0.5 )
38 changes: 23 additions & 15 deletions BernGrid.R
View file
Original file line number Diff line number Diff line change
Expand Up @@ -25,7 +25,8 @@ BernGrid = function( Theta , pTheta , Data ,
# Hints:
# You will need to "source" this function before calling it.
# You may want to define a tall narrow window before using it; e.g.,
# > windows(7,10)
# > source("openGraphSaveGraph.R")
# > openGraph(width=7,height=10,mag=0.7)

# Create summary values of Data
z = sum( Data==1 ) # number of 1's in Data
Expand All @@ -39,7 +40,7 @@ pThetaGivenData = pDataGivenTheta * pTheta / pData
# Plot the results.
layout( matrix( c( 1,2,3 ) ,nrow=3 ,ncol=1 ,byrow=FALSE ) ) # 3x1 panels
par( mar=c(3,3,1,0) , mgp=c(2,1,0) , mai=c(0.5,0.5,0.3,0.1) ) # margin settings
dotsize = 4 # how big to make the plotted dots
dotsize = 5 # how big to make the plotted dots
# If the comb has a zillion teeth, it's too many to plot, so plot only a
# thinned out subset of the teeth.
nteeth = length(Theta)
Expand All @@ -54,7 +55,7 @@ meanTheta = sum( Theta * pTheta ) # mean of prior, for plotting
plot( Theta[thinIdx] , pTheta[thinIdx] , type="p" , pch="." , cex=dotsize ,
xlim=c(0,1) , ylim=c(0,1.1*max(pThetaGivenData)) , cex.axis=1.2 ,
xlab=bquote(theta) , ylab=bquote(p(theta)) , cex.lab=1.5 ,
main="Prior" , cex.main=1.5 )
main="Prior" , cex.main=1.5 , col="skyblue" )
if ( meanTheta > .5 ) {
textx = 0 ; textadj = c(0,1)
} else {
Expand All @@ -68,7 +69,7 @@ plot(Theta[thinIdx] ,pDataGivenTheta[thinIdx] ,type="p" ,pch="." ,cex=dotsize
,xlim=c(0,1) ,cex.axis=1.2 ,xlab=bquote(theta)
,ylim=c(0,1.1*max(pDataGivenTheta))
,ylab=bquote( "p(D|" * theta * ")" )
,cex.lab=1.5 ,main="Likelihood" ,cex.main=1.5 )
,cex.lab=1.5 ,main="Likelihood" ,cex.main=1.5 , col="skyblue" )
if ( z > .5*N ) { textx = 0 ; textadj = c(0,1) }
else { textx = 1 ; textadj = c(1,1) }
text( textx ,1.0*max(pDataGivenTheta) ,cex=2.0
Expand All @@ -78,7 +79,7 @@ meanThetaGivenData = sum( Theta * pThetaGivenData )
plot(Theta[thinIdx] ,pThetaGivenData[thinIdx] ,type="p" ,pch="." ,cex=dotsize
,xlim=c(0,1) ,ylim=c(0,1.1*max(pThetaGivenData)) ,cex.axis=1.2
,xlab=bquote(theta) ,ylab=bquote( "p(" * theta * "|D)" )
,cex.lab=1.5 ,main="Posterior" ,cex.main=1.5 )
,cex.lab=1.5 ,main="Posterior" ,cex.main=1.5 , col="skyblue" )
if ( meanThetaGivenData > .5 ) { textx = 0 ; textadj = c(0,1) }
else { textx = 1 ; textadj = c(1,1) }
text(textx ,1.00*max(pThetaGivenData) ,cex=2.0
Expand All @@ -94,16 +95,23 @@ points( Theta[ HDIinfo$indices ] ,
text( mean( Theta[ HDIinfo$indices ] ) , HDIinfo$height ,
bquote( .(100*signif(HDIinfo$mass,3)) * "% HDI" ) ,
adj=c(0.5,-1.5) , cex=1.5 )
# Mark the left and right ends of the waterline. This does not mark
# internal divisions of an HDI waterline for multi-modal distributions.
lowLim = Theta[ min( HDIinfo$indices ) ]
highLim = Theta[ max( HDIinfo$indices ) ]
lines( c(lowLim,lowLim) , c(-0.5,HDIinfo$height) , type="l" , lty=2 , lwd=1.5)
lines( c(highLim,highLim) , c(-0.5,HDIinfo$height) , type="l" , lty=2 , lwd=1.5)
text( lowLim , HDIinfo$height , bquote(.(round(lowLim,3))) ,
adj=c(1.1,-0.1) , cex=1.2 )
text( highLim , HDIinfo$height , bquote(.(round(highLim,3))) ,
adj=c(-0.1,-0.1) , cex=1.2 )
# Mark the left and right ends of the waterline.
# Find indices at ends of sub-intervals:
inLim = HDIinfo$indices[1] # first point
for ( idx in 2:(length(HDIinfo$indices)-1) ) {
if ( ( HDIinfo$indices[idx] != HDIinfo$indices[idx-1]+1 ) | # jumps on left, OR
( HDIinfo$indices[idx] != HDIinfo$indices[idx+1]-1 ) ) { # jumps on right
inLim = c(inLim,HDIinfo$indices[idx]) # include idx
}
}
inLim = c(inLim,HDIinfo$indices[length(HDIinfo$indices)]) # last point
# Mark vertical lines at ends of sub-intervals:
for ( idx in inLim ) {
lines( c(Theta[idx],Theta[idx]) , c(-0.5,HDIinfo$height) , type="l" , lty=2 ,
lwd=1.5 )
text( Theta[idx] , HDIinfo$height , bquote(.(round(Theta[idx],3))) ,
adj=c(0.5,-0.1) , cex=1.2 )
}

return( pThetaGivenData )
} # end of function
31 changes: 31 additions & 0 deletions BernGridExample.R
View file
Original file line number Diff line number Diff line change
@@ -0,0 +1,31 @@
graphics.off()
source("openGraphSaveGraph.R") # for openGraph() function, used below.
source("BernGrid.R")

# For Figure 6.4:
# Specify theta values.
thetagrid = seq(0,1,length=1001)
# Specify probability mass at each theta value.
relprob = sin( 2*pi*thetagrid )^6
prior = relprob / sum(relprob) # probability mass at each theta
# Specify the data vector.
datavec = c( rep(1,2) , rep(0,1) )
# Open a window.
openGraph(width=7,height=10,mag=0.7)
# Call the function.
posterior = BernGrid( Theta=thetagrid , pTheta=prior , Data=datavec )
saveGraph(file="Fig.6.4",type="jpg")

# For Figure 6.5:
pTheta = c( 50:1 , rep(1,50) , 1:50 , 50:1 , rep(1,50) , 1:50 )
pTheta = pTheta / sum( pTheta )
width = 1 / length(pTheta)
Theta = seq( from = width/2 , to = 1-width/2 , by = width )
dataVec = c( rep(1,3) , rep(0,1) )
openGraph(width=7,height=10,mag=0.7)
posterior = BernGrid( Theta=Theta , pTheta=pTheta , Data=dataVec )
saveGraph(file="Fig.6.5left",type="jpg")
dataVec = c( rep(1,12) , rep(0,4) )
openGraph(width=7,height=10,mag=0.7)
posterior = BernGrid( Theta=Theta , pTheta=posterior , Data=dataVec )
saveGraph(file="Fig.6.5right",type="jpg")
Loading

0 comments on commit fa4c5ba

Please sign in to comment.