# Normal approximation to posterior distribution
# (aka, the Bayesian Central Limit Theorem)
# Discussed in class on Tue Feb 23
rm(list=ls())
par(cex.lab=1.5,cex.main=1.5,cex.lab=1.5,mar=c(5,4+0.5,4,2)+0.1)
# Small n case
theta <- seq(0.001,0.999,.001) # Grid of thetas for plotting n <- 10; y <- 3; # The data # Compute the posterior mean and Fisher information matrix A <- y B <- n - y theta.MAP <- A / (A+B) Info <- A/theta.MAP^2 + B/(1-theta.MAP)^2 # Plot the true and approximate posteriors post1 <- dbinom(y,n,theta) * dbeta(theta,1,1) post1 <- post1/sum(post1) post2 <- dnorm(theta, theta.MAP, sqrt(1/Info)) post2 <- post2/sum(post2) plot(theta, post1, type="l", lwd=2, xlab=expression(theta), ylab="Posterior", main="n=10, y=3") abline(v=theta.MAP, lty=3, lwd=2) lines(theta, post2, lty=2, lwd=2) legend("topright",c("Exact","CLT","MAP"),bty="n",lty=1:3,cex=1.5,lwd=2) # Medium n case rm(list=ls()) theta <- seq(0.001,0.999,.001) # Grid of thetas for plotting n <- 30; y <- 9; # The data # Compute the posterior mean and Fisher information matrix A <- y B <- n - y theta.MAP <- A/(A+B) Info <- A/theta.MAP^2 + B/(1-theta.MAP)^2 # Plot the true and approximate posteriors post1 <- dbinom(y,n,theta) * dbeta(theta,1,1) post1 <- post1/sum(post1) post2 <- dnorm(theta, theta.MAP, sqrt(1/Info)) post2 <- post2/sum(post2) plot(theta, post1, type="l", lwd=2, xlab=expression(theta),ylab="Posterior", main="n=30, y=9") abline(v=theta.MAP,lty=3,lwd=2) lines(theta, post2, lty=2, lwd=2) legend("topright",c("Exact","CLT","MAP"),bty="n",lty=1:3,cex=1.5,lwd=2) # Approximation in a large sample case rm(list=ls()) theta <- seq(0.001,0.999,.001) # Grid of thetas for plotting n <- 100; y <- 30; # The data # Compute the posterior mean and Fisher information matrix A <- y B <- n - y theta.MAP <- A/(A+B) Info <- A/theta.MAP^2 + B/(1-theta.MAP)^2 # Plot the true and approximate posteriors post1 <- dbinom(y,n,theta) * dbeta(theta,1,1) post1 <- post1/sum(post1) post2 <- dnorm(theta, theta.MAP, sqrt(1/Info)) post2 <- post2/sum(post2) plot(theta, post1, type="l", lwd=2, xlab=expression(theta), ylab="Posterior", main="n=100, y=30") abline(v=theta.MAP, lty=3, lwd=2) lines(theta, post2, lty=2, lwd=2) legend("topright",c("Exact","CLT","MAP"),bty="n",lty=1:3,cex=1.5,lwd=2)