######################################
# Simulated data
######################################
set.seed(1234)
n  = 10
x = rnorm(n,0,1)
#x = rt(n,4)/sqrt(2)
x

######################################
# model 0
# x1,...,xn iid N(theta,1)
######################################

# posterior mean and posterior variance
mtheta = sum(x)/(n+1)
vtheta = 1/(n+1)
c(mtheta,sqrt(vtheta))

# posterior 95% credibility interval
qnorm(c(0.025,0.975),mtheta,sqrt(vtheta))

# prior predictive
pred0=(2*pi)^(-n/2)*sqrt(1/(1+n))*exp(-0.5*sum(x^2))*exp(0.5/(1+n)*(sum(x))^2)

pred0

######################################
# model 0
# x1,...,xn iid N(theta,1)
######################################

# prior predictive (via Monte Carlo integration)
set.seed(1234)
A = gamma(2.5)/gamma(2)/sqrt(2*pi)
B = A^n/sqrt(2*pi)

M = 1000000
thetas = rnorm(M,0,sqrt(10))
hist(thetas)

g = function(theta){
  exp(-theta^2/2)*prod((1+(x-theta)^2/2)^(-2.5))
}

int = 0
for (j in 1:M)
  int = int + g(thetas[j])/dnorm(thetas[j],0,sqrt(10))
pred=int/M
pred1 = B*pred

# Bayes factor of model 0 versus model 1
c(pred0,pred1,pred0/pred1)


# Comparing posterior densities of theta
thetas1 = seq(-1.5,1,length=1000)
post1 = rep(0,1000)
for (j in 1:1000)
  post1[j]= B*g(thetas1[j])/pred1

plot(thetas1,dnorm(thetas1,mtheta,sqrt(vtheta)),
type="l",ylim=c(0,1.5),xlab=expression(theta),ylab="Posterior density")
lines(thetas1,post1,col=2)
legend("topleft",legend=c("Model 0","Model 1"),col=1:2,lty=1)