Code for producing sequential updating plots for globe tossing example in Chapter 2 of Statistical Rethinking (Figure 2.5, page 30)

ch2_globe_toss_bayes_update.R

# chapter 2 code

library(rethinking)

col1 <- "black"
col2 <- col.alpha( "black" , 0.7 )

#####
# introducing the golem
# show posterior as data comes in

# 1 indicates 'water'; 0 indicates 'land'
d <- c(1,0,1,1,1,0,1,0,1) # length 9
l <- ifelse( d==1 , "W" , "L" )

# prior uniform - "Laplace" prior
a <- 1
b <- 1

# draw first plot
i <- 0
curve( dbeta(x,a,b) , from=0 , to=1 , xlab="proportion water" , ylab="" , col=col1 , lwd=2 , ylim=c(0,2.75) , xaxt="n" , yaxt="n" )
axis( 1 , at=c(0,0.5,1) , labels=c(0,0.5,1) )
title( ylab="plausibility" , mgp=c(1,1,0) )
text( 0.15 , 2.5 , paste("n =",i) )
mtext( paste(l) , at=seq(from=0,to=1,length.out=length(l)) , col=c( rep( "darkgray" ,length(l)) ) )

# execute this code once for each subsequent plot
i <- i + 1
prior.a <- a
prior.b <- b
a <- a + d[i]
b <- b + 1 - d[i]
curve( dbeta(x,a,b) , from=0 , to=1 , xlab="probability of water" , ylab="" , col=col1 , lwd=2 , ylim=c(0,2.75) , xaxt="n" , yaxt="n" )
axis( 1 , at=c(0,0.5,1) , labels=c(0,0.5,1) )
text( 0.15 , 2.5 , paste("n =",i) )
curve( dbeta(x,prior.a,prior.b) , add=TRUE , col=col1 , lty=2 )
mtext( paste(l) , at=seq(from=0,to=1,length.out=length(l)) , col=c( rep(col1,i) , rep( "darkgray" ,length(l)-i) ) )