70 CH A P T E R 4: Bayes' Rule The listing that follows includes line numbers in the margins, to facilitate track- ing the code across page splits and to facilitate referring to specific lines of the code when you have enthusiastic conversations about it at parties. Mac users: If you are running R under MacOS instead of in a Windows emu- lator such as WINE, you will need to change all the windows() commands to quartz() . Later in the book, when we use BUGS, there is no Mac equivalent and you must run the programs under WINE or windows. (BayesUpdate.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 # Theta is the vector of candidate values for the parameter theta. # nThetaVals is the number of candidate theta values. # To produce the examples in the book, set nThetaVals to either 3 or 63. nThetaVals = 3 # Now make the vector of theta values: Theta = seq( from = 1/(nThetaVals+1) , to = nThetaVals/(nThetaVals+1) , by = 1/(nThetaVals+1) ) } # pTheta is the vector of prior probabilities on the theta values. 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 ) # Compute the likelihood of the data for each value of theta: pDataGivenTheta = Theta^nHeads * (1-Theta)^nTails # Compute the posterior: 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. 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 par(mai = c(0.5,0.5,0.3,0.1)) # margin size in inches: bottom,left,top,right # Plot the prior: 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 ) # Plot the likelihood: plot( Theta , pDataGivenTheta , type = "h" , lwd = 3 , main = "Likelihood" , xlim = c(0,1) , xlab = bquote(theta) ,