Model

We will illustrate the use of Monte Carlo integration via a simple normal-normal example where observations are iid normal, ie. \[ y_1,\ldots,y_n|\theta \sim N(\theta,1) \] and the prior information regarding \(\theta\) is \[ \theta \sim N(0,1) \]

Simulating the dataset

set.seed(23235)
n = 10 
thetatrue = 1
y = rnorm(n,thetatrue,1)
ybar = mean(y)
ybar
## [1] 1.350405

In this simple case, all computations can be done in analytical form. More precisely, it is easy to see that the sampling distribution of the sample mean \({\bar y}_n=\sum_{i=1}^n y_i\), given \(\theta\), is \[ {\bar y}_n|\theta \sim N(\theta,1/n) \] while the prior predictive of \({\bar y}_n\) is \[ {\bar y}_n \sim N\left(0,\frac{n+1}{n}\right) \] Similarly, the posterior distribution of \(\theta\) given \({\bar y}_n\) is also Gaussian: \[ \theta | {\bar y}_n \sim N\left(\frac{n}{n+1}{\bar y}_n,\frac{1}{n+1}\right) \] and the posterior predictive of \(y_{n+1}\) is \[ p(y_{n+1}|{\bar y}_n) \sim N\left(\frac{n}{n+1}{\bar y}_n,\frac{n+2}{n+1}\right) \]

Prior predictive of \({\bar y}_n\)

par(mfrow=c(1,1))
yy = seq(-5,5,length=100)
plot(yy,dnorm(yy,0,sqrt((n+1)/n)),type="l",ylab="Density",xlab=expression(bar(y)))
legend("topleft",legend=c("Model",expression(y[i]~"|"~theta %~%N(theta,1)~i==list(1,...,n)),expression(bar(y)%~%N(0,1/n))))
legend("topright",legend=c("Prior",expression(theta %~%N(0,1))))
title(substitute(paste("Prior predictive of ",bar(y)," (n = ",n,")",sep=""),list(n=n)))
points(ybar,0,col=2,pch=16,cex=1)

Posterior of \(\theta\)

m = n/(n+1)*ybar
sd = sqrt(1/(n+1))
tt = seq(m-5*sd,m+5*sd,length=100)
plot(tt,dnorm(tt,m,sd),xlab=expression(theta),ylab="Density",type="l")
points(thetatrue,0,col=2,pch=16,cex=1)
title(expression('Posterior of ' ~ theta))

Prior and posterior predictives

m = n/(n+1)*ybar
sd = sqrt((n+2)/(n+1))
yy1 = seq(m-10*sd,m+10*sd,length=1000)
plot(yy1,dnorm(yy1,m,sd),xlab=expression(y[n+1]),ylab="Density",type="l")
lines(yy1,dnorm(yy1,0,2),col=2)
legend("topleft",legend=c(expression(Pr(y[n+1])),expression(Pr(y[n+1]~"|"~list(y[1],...,y[n])))),col=2:1,lty=1)
legend("topright",legend=substitute(paste(bar(y)==yyy),list(yyy=round(ybar,3))))