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)
{
## Simulate outcomes in advance
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)
 curve(dbeta(x,b_a,b_b),col='blue',lty=2,add=TRUE)

 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)
 #{
 # curve(dbeta(x,a_a+success[i]+1,a_b+(i-success[i])+1),add=TRUE,col=rgb(0,100,0,(1-(frame-i)/frame) * 100,maxColorValue=255))
 # curve(dbeta(x,b_a+success[i]+1,b_b+(i-success[i])+1),add=TRUE,col=rgb(0,0,100,(1-(frame-i)/frame) * 100,maxColorValue=255))
 #}

 curve(dbeta(x,a_a+success[i]+1,a_b+(i-success[i])+1),add=TRUE,col=rgb(0,100,0,255,maxColorValue=255),lwd=2)

 curve(dbeta(x,b_a+success[i]+1,b_b+(i-success[i])+1),add=TRUE,col=rgb(0,0,100,255,maxColorValue=255),lwd=2)

 # 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
# http://www.youtube.com/watch?v=u5_3MGP2Lj4&amp;feature=player_detailpage#t=265s
}
sim_bayes(a_a=1,a_b=1,b_a=20,b_b=20,p=0.25,N=250)
Advertisements

2 thoughts on “My take on the Bayesian updating video

    1. Sam Clifford Post author

      Thanks, Corey! I ran through my example with one of my coworkers this morning. He’s not a statistician let alone a Bayesian but I did manage to use the example to explain how Bayesian updating works and what sort of influence the prior can have on the results.

      Reply

Leave a Reply

Fill in your details below or click an icon to log in:

WordPress.com Logo

You are commenting using your WordPress.com account. Log Out / Change )

Twitter picture

You are commenting using your Twitter account. Log Out / Change )

Facebook photo

You are commenting using your Facebook account. Log Out / Change )

Google+ photo

You are commenting using your Google+ account. Log Out / Change )

Connecting to %s