############################################################
# A) Comparing MCI and MC via importance function
#
# E(theta) where theta ~ Beta(8,4)
# Exact computation: E(theta) = 8/12 = 2/3
############################################################

# Monte Carlo integration
M   = 10000
sim = matrix(0,20,M)
for (i in 1:20){
  theta = rbeta(M,8,4)
  gtheta = theta
  sim[i,] = cumsum(gtheta)/(1:M)
}

# Monte Carlo integration via importance function
# Importance function: U(a,b)
a    = 0.3
b    = 1.0
sim1 = matrix(0,20,M)
for (i in 1:20){
  theta = runif(M,a,b)
  gtheta = theta*dbeta(theta,8,4)/(1/(b-a))
  sim1[i,] = cumsum(gtheta)/(1:M)
}


par(mfrow=c(1,1))
ts.plot(t(sim1),col=grey(0.75),ylim=c(0.6,0.75),
xlab="Monte Carlo sample size",ylab="Expectation")
for (i in 1:20)
  lines(sim[i,])
legend("topright",legend=c("Exact integral","Monte Carlo Integration","MC via importance funtion"),col=c(2,grey(0.75),1),lwd=2,bty="n")
abline(h=2/3,col=2,lwd=2)


############################################################
# B) Raoblackwellization in action
#
# Approximating E(Y) via MCI with draws from Y
# Approximating E(E(Y|X)) via MCI with draws from X
#
############################################################
# library("mvtnorm")
M   = 1000
mu  = c(0,1)
Sig = matrix(c(1,0.25,0.25,1),2,2)
R   = 1000
I1  = rep(0,R)
I2  = rep(0,R)
for (r in 1:R){
  draws = rmvnorm(M,mu,Sig)
  I1[r] = mean(draws[,2])
  I2[r] = mu[2]+Sig[1,2]/Sig[1,1]*(mean(draws[,1])-mu[1])
}
boxplot(I1,I2,outline=FALSE)