My take on the Bayesian updating video

I’m referring here to my last post where I reblogged something that Bayesian Biologist had written. I’ve modified the code (and included it below) to tweak it to my tastes and show the different behaviour of the posterior under a uniform Beta(1,1) prior and an informative Beta(20,20) prior. Rather than flipping a fair coin I’ve assumed that the coin is unfair (p=0.25) and the Beta(20,20) prior is chosen by the naive observer who believes, with a high degree of certainty, that the coin is fair.

This is, to me, an illustration of the power of uninformative priors (something Laplace seemed keen on in Binomial experiments). It’s also a good demonstration of how given enough data the posteriors will converge to the same result. A persistent criticism of Bayesian analysis is that the priors are so subjective that anyone can come up with their own prior and get a different result to another observer. As more data is collected, the influence of the prior in the posterior is diminished.

We see (in the 1.08MB animation below the cut) that the Beta(1,1) prior is more sensitive to each success and failure and the posterior mean approximates p much quicker than the more “certain” Beta(20,20) prior. The take home message? Be more vague than you think is warranted just in case your prior is not diffuse enough.

To generate the gif I’ve used the instructions in the video in the last lines of the code.

## Corey Chivers, 2012 ##
## modifications by ##
## Sam Clifford, 2012 ##
sim_bayes <- function(p=0.5,N=100,y_lim=20,a_a=2,a_b=10,b_a=8,b_b=3)
{
outcomes<-sample(1:0,N,prob=c(p,1-p),replace=TRUE)
success<-cumsum(outcomes)

for(frame in 1:N)
{
png(paste("plots/",1000+frame,".png",sep=""))
curve(dbeta(x,a_a,a_b),xlim=c(0,1),ylim=c(0,y_lim),col='green',xlab='p',ylab='Posterior Density',lty=2)

lines(x=c(p,p),y=c(0,y_lim),lty=2,lwd=1,col="grey60")

i <- frame
# i don't like having the leftovers on the screen -- Sam
#for(i in 1:frame)
#{
#}

# modifications to make prior explicit in legend
legend('topleft',legend=c( paste('Observer A - Beta(', a_a , ',' , a_b , ')',sep=""),paste('Observer B - Beta(', b_a , ',' , b_b , ')',sep="")),lty=1,col=c('green','blue'))
text(0.75,17,label=paste(success[i],"successes,",i-success[i],"failures"))

dev.off()
}
# Sam
# system('mencoder mf://plots/*.png -mf fps=15:type=png -ovc copy -oac copy -o plots/output.avi')
# we'll use GIMP rather than mencoder