Turning the Bayesian crank

Normal data and four priors

1 Model and priors

In this example we use one unique model structure \[ x|\theta \sim N(\theta,\sigma^2) \] where \(\sigma=40\). The goal is to learn about \(\theta\) from a Bayesian view point where four alternative prior distributions are considered: \[\begin{eqnarray*} \mbox{Prior A} &:& \theta \sim N(\mu_A,\sigma_A^2)\\ \mbox{Prior B} &:& \theta \sim N(\mu_B,\sigma_B^2)\\ \mbox{Prior C} &:& \theta \sim t_\nu(\mu_C,\sigma_C^2)\\ \mbox{Prior D} &:& \theta \sim \frac{1}{3}N(\mu_A,\sigma_A^2)+ \frac{1}{3}N(\mu_A,\sigma_A^2)+ \frac{1}{3}t_\nu(\mu_C,\sigma_C^2), \end{eqnarray*}\] and prior hyperparameters: \[\begin{eqnarray*} \mu_A &=& 900,\sigma_A=20\\ \mu_B &=& 800,\sigma_B=80\\ \mu_C &=& 900,\sigma_C=\sqrt{240}, \nu=5. \end{eqnarray*}\]

den.t = function(theta,mu,sig,df){
  dt((theta-mu)/sig,df=df)/sig
}
prior.D = function(theta){
  (dnorm(theta,mu.A,sig.A)+
  dnorm(theta,mu.B,sig.B)+
  den.t(theta,mu.C,sig.C,nu))/3
}

sig = 40
x = 850

mu.A= 900
sig.A = 20

mu.B = 800
sig.B = 80

mu.C = 900
sig.C = sqrt(240)
nu     = 5

2 Prior densities

par(mfrow=c(1,1))
theta = seq(500,1100,length=1000)
plot(theta,dnorm(theta,x,sig),type="l",lwd=2,ylim=c(0,0.025),
        xlab=expression(theta),ylab="Priors & likelihood")
lines(theta,dnorm(theta,mu.A,sig.A),col=2,lwd=2)
lines(theta,dnorm(theta,mu.B,sig.B),col=3,lwd=2)
lines(theta,den.t(theta,mu.C,sig.C,nu),col=4,lwd=2)
lines(theta,prior.D(theta),col=5,lwd=2)
points(x,0,pch=16,cex=2)
legend("topleft",legend=c("Likelihood","Prior A","Prior B",
       "Prior C","Prior D"),col=1:5,lwd=2)

3 Logarithm of the prior densities

par(mfrow=c(1,1))
plot(theta,dnorm(theta,x,sig,log=TRUE),type="l",lwd=2,ylim=c(-20,0),
        xlab=expression(theta),ylab="Log priors & log likelihood")
lines(theta,dnorm(theta,mu.A,sig.A,log=TRUE),col=2,lwd=2)
lines(theta,dnorm(theta,mu.B,sig.B,log=TRUE),col=3,lwd=2)
lines(theta,log(den.t(theta,mu.C,sig.C,nu)),col=4,lwd=2)
lines(theta,log(prior.D(theta)),col=5,lwd=2)
legend("topleft",legend=c("Likelihood","Prior A","Prior B",
       "Prior C","Prior D"),col=1:5,lwd=2,bty="n")

4 Posterior densities via SIR

Posterior summary via sampling importance resampling (SIR)

set.seed(125532)
M      = 1000000
thetas = runif(M,600,1000)
w      = dnorm(thetas,x,sig)*dnorm(thetas,mu.A,sig.A)
ind.A  = sample(1:M,size=M,replace=TRUE,prob=w)
w      = dnorm(thetas,x,sig)*dnorm(thetas,mu.B,sig.B)
ind.B  = sample(1:M,size=M,replace=TRUE,prob=w)
w      = dnorm(thetas,x,sig)*den.t(thetas,mu.C,sig.C,nu)
ind.C  = sample(1:M,size=M,replace=TRUE,prob=w)
w      = dnorm(thetas,x,sig)*(dnorm(thetas,mu.A,sig.A)+
         dnorm(thetas,mu.B,sig.B)+den.t(thetas,mu.C,sig.C,nu))
ind.D  = sample(1:M,size=M,replace=TRUE,prob=w)

par(mfrow=c(1,1))
plot(density(thetas[ind.A]),type="l",xlab=expression(theta),
        ylab="Density",ylim=c(0,0.025),xlim=c(650,950),lwd=2,main="")
lines(density(thetas[ind.B]),col=2,lwd=2)
lines(density(thetas[ind.C]),col=4,lwd=2)
lines(density(thetas[ind.D]),col=3,lwd=2)
points(x,0,pch=16,cex=2,col=6)
legend("topleft",legend=c("Posterior A","Posterior B","Posterior C","Posterior D"),col=c(1,2,4,3),lwd=2)
title(paste("Data point x = ",x,sep=""))

5 Posterior summaries

Below we us Monte Carlo approximations to approximate \[\begin{eqnarray*} E(\theta|x)\\ Median(\theta|x)\\ Pr(\theta < x|x)\\ Pr(\theta<q_{0.05})=0.05. \end{eqnarray*}\]

tab = rbind(

c(mean(thetas[ind.A]),mean(thetas[ind.B]),
mean(thetas[ind.C]),mean(thetas[ind.D])),

c(median(thetas[ind.A]),median(thetas[ind.B]),
median(thetas[ind.C]),median(thetas[ind.D])),

c(mean(thetas[ind.A]<x),mean(thetas[ind.B]<x),
mean(thetas[ind.C]<x),mean(thetas[ind.D]<x)),

c(quantile(thetas[ind.A],0.05),quantile(thetas[ind.B],0.05),
quantile(thetas[ind.C],0.05),quantile(thetas[ind.D],0.05)))

rownames(tab) = c("E(theta)","Median(theta)",
                  "Pr(theta<x)","Quantile(0.05)")

colnames(tab) = c("Post A","Post B","Post C","Post D")
round(tab,2)

##                Post A Post B Post C Post D
## E(theta)       890.06 840.02 891.16 876.15
## Median(theta)  890.10 840.09 892.45 884.35
## Pr(theta<x)      0.01   0.61   0.02   0.19
## Quantile(0.05) 860.52 781.18 860.30 806.53