##################################################################
#
# Monte Carlo integration - 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/
#
#################################################################

# Computing the integral from 0 to 1 of 
# h(theta) = (cos(50*theta)+sin(20*theta))^2
# via simple Monte Carlo with draws from U(0,1)
Ms    = c(100,1000,10000,100000,1000000)
hbars = NULL
ints  = NULL
for (M in Ms){
  theta = runif(M,0,1)
  h     = (cos(50*theta)+sin(20*theta))^2
  hbar  = sum(h)/M
  var   = sum((h-hbar)^2)/M
  mc.se = sqrt(var/M)
  int   = c(hbar-2*mc.se,hbar+2*mc.se)
  hbars = c(hbars,hbar)
  ints  = rbind(ints,int)
}
cbind(Ms,hbars,ints)

# Computing Pr(theta>2) where theta ~ Cauchy(0,1)
# Simple Monte Carlo: sampling from Cauchy(0,1)
Ms    = c(100,1000,10000,100000,1000000)
hbars = NULL
ints  = NULL
for (M in Ms){
  theta = rt(M,1)
  h = rep(0,M)
  h[theta>2]=1
  hbar = sum(h)/M
  var = sum((h-hbar)^2)/M
  mc.se = sqrt(var/M)
  int = c(hbar-2*mc.se,hbar+2*mc.se)
  hbars = c(hbars,hbar)
  ints = rbind(ints,int)
}
simplemc = cbind(Ms,hbars,ints)

# Computing Pr(theta>2) where theta ~ Cauchy(0,1)
# by using the property: Pr(theta>2)=0.5*Pr(|theta|<2)
# Simple Monte Carlo: sampling from Cauchy(0,1)
Ms    = c(100,1000,10000,100000,1000000)
hbars = NULL
ints  = NULL
for (M in Ms){
  theta = rt(M,1)
  h = rep(0,M)
  h[abs(theta)>2]=1
  hbar = sum(h)/M/2
  var = sum((h-hbar)^2)/M
  mc.se = sqrt(var/M)
  int = c(hbar-2*mc.se,hbar+2*mc.se)
  hbars = c(hbars,hbar)
  ints = rbind(ints,int)
}
notsosimplemc = cbind(Ms,hbars,ints)
  
# Computing Pr(theta>2) where theta ~ Cauchy(0,1)
# by using the property: Pr(theta>2) = integral from 
# 0 to 1/2 of (1/u^2)/(2*pi*(1+1/u^2))
# Simple Monte Carlo: sampling from U(0,1/2)
Ms    = c(100,1000,10000,100000,1000000)
hbars = NULL
ints  = NULL
for (M in Ms){
  theta = runif(M,0,1/2)
  h    = (1/theta^2)/(2*pi*(1+1/theta^2))
  hbar = sum(h)/M
  var = sum((h-hbar)^2)/M
  mc.se = sqrt(var/M)
  int = c(hbar-2*mc.se,hbar+2*mc.se)
  hbars = c(hbars,hbar)
  ints = rbind(ints,int)
}
prettyclevermc = cbind(Ms,hbars,ints)


# Rejection method for the normal normal example
pi = function(theta){
  exp(-0.5*(x-theta)^2)*exp(-0.5*(theta-theta0)^2/C0)
}
q = function(theta){
  dt((theta-x)/5,10)
}

x      = 1
theta0 = 0
C0     = 1

# What is A?
thetas = seq(-30,30,length=1000)
par(mfrow=c(1,1))
plot(thetas,ps,type="l",ylim=c(0,0.85))
lines(thetas,2.1*qs,col=2,lwd=3)

# Drawing from proposal q(.)
thetas  = x+sqrt(5)*rt(M,10)

# Computing acceptance probabilities
prob    = pi(thetas)/(2.1*q(thetas))

# Drawing the uniforms
u       = runif(M)

# Draws from pi(.) (the target distribution)
thetas1 = thetas[u<prob] 

par(mfrow=c(1,1))
hist(thetas1,prob=TRUE,main="")
ths = seq(x/2-3*sqrt(1/2),x/2+3*sqrt(1/2),length=1000)
lines(ths,dnorm(ths,x/2,sqrt(1/2)),col=2,lwd=2)
title(paste("Acceptance=",length(thetas1)/M,sep=""))


# Sampling importance resampling

# What I want is to draw from N(0,1)
M = 10000
theta = rnorm(M,0,2)

w = dnorm(theta,0,1)/dnorm(theta,0,2)


plot(theta,w,type="h")

theta1 = sample(theta,replace=TRUE,size=M,prob=w)

hist(theta1,prob=TRUE)
grid = seq(min(theta1),max(theta1),length=1000)
lines(grid,dnorm(grid,0,1),col=2,lwd=3)
length(unique(theta1))

# Previous normal-normal example using SIR
x      = 1
theta0 = 0
C0     = 1
M      = 10000
theta  = rnorm(M,theta0,sqrt(C0))
w      = dnorm(x,theta,1)
theta1 = sample(theta,replace=TRUE,size=M,prob=w)

hist(theta1,prob=TRUE)
grid = seq(min(theta1),max(theta1),length=1000)
lines(grid,dnorm(grid,x/2,sqrt(1/2)),col=2,lwd=3)
length(unique(theta1))