##################################################################
# SIR - classroom
#
##################################################################
# HEDIBERT FREITAS LOPES
# Associate Professor of Econometrics and Statistics
# The University of Chicago Booth School of Business
# 5807 South Woodlawn Avenue
# Chicago, Illinois, 60637
# Email : hlopes@ChicagoGSB.edu
# URL: http://faculty.chicagobooth.edu/hedibert.lopes/research/
#
#################################################################
q1 = function(theta){
  dt((theta-x)/sqrt(V1),nu)
}
q2 = function(theta){
  dnorm(theta,theta0,sqrt(C0))
}
post = function(theta){
  dnorm(theta,E.true,sqrt(V.true))
}
like = function(theta){
  dnorm(x,theta,1)
}
prior = function(theta){
  dnorm(theta,theta0,sqrt(C0))
}

nu     = 10
theta0 = 0
C0     = 1
V1     = 5
x      = 1
E.true = (theta0+x*C0)/(1+C0)
V.true = C0/(1+C0)


ES  = NULL
ESr = NULL
M = 1000
for (i in 1:100){
  # Sampling from q1 and q2
  theta.q1 =      x + sqrt(V1)*rt(M,nu)
  theta.q2 = theta0 + sqrt(C0)*rnorm(M)

  w1 = like(theta.q1)*prior(theta.q1)/q1(theta.q1)
  w2 = like(theta.q2)*prior(theta.q2)/q2(theta.q2)

  E1 = sum(w1*theta.q1)/sum(w1)
  E2 = sum(w2*theta.q2)/sum(w2)
  
  theta1r = sample(theta.q1,size=M,replace=TRUE,prob=w1)
  theta2r = sample(theta.q2,size=M,replace=TRUE,prob=w2)
  
  E1r = mean(theta1r)
  E2r = mean(theta2r)
  
  ES  = rbind(ES,c(E1,E2))
  ESr = rbind(ESr,c(E1r,E2r))
  
  
}

L = min(ES,ESr)
U = max(ES,ESr)
par(mfrow=c(1,2))
boxplot(ES[,1],ES[,2],names=c("Student's t","Prior"),ylim=c(L,U))
abline(h=E.true,col=3,lwd=4)
title("Reweighting")

boxplot(ESr[,1],ESr[,2],names=c("Student's t","Prior"),ylim=c(L,U))
abline(h=E.true,col=3,lwd=4)
title("Resampling")