Monte Carlo Integration and Simulation

Five illustrative examples

PROBLEM 1: Grid approximation

Here we wanto to compute \(E(\theta)\), \(E(\theta^2)\) and \(Pr(\theta<0)\), when the unnormalized density of \(\theta\) is \[ p(\theta) \propto \exp\{-0.5(\theta-1)^2\}(1+0.25\theta^2)^{-2.5} \qquad \theta \in \Re. \] We construct a grid from \(-5\) to \(5\) with one million points, so the distance between two consecutive points is about \(0.00001\). The figure below show that the variation of \(\theta\) is roughly between \(-4\) and \(4\), so our grid is quite conservative. We compare the performance of the grid approximation for various sizes of the grid, from 10 points to 1,000,000 points.

theta = seq(-5,5,length=1000000)
h     = theta[2]-theta[1]
post  = exp(-0.5*(theta-1)^2)*(1+0.25*theta^2)^(-2.5)
plot(theta,post,type="l",xlab=expression(theta),
     ylab="Unnormalized density")

Ms = c(seq(10,100,by=10),1000,10000,100000,1000000)
summaries = NULL
for (M in Ms){
  theta = seq(-5,5,length=M)
  post  = exp(-0.5*(theta-1)^2)*(1+0.25*theta^2)^(-2.5)
  denominator = sum(post)
  E = sum(theta*post)/denominator
  E2 = sum((theta-E)^2*post)/denominator
  Prob = sum(ifelse(theta<0,1,0)*post)/denominator
  summaries=rbind(summaries,round(c(E,E2,Prob),10))
}
tab = cbind(Ms,summaries)
colnames(tab) = c("M","Mean","Mean2","Prob(theta<0)")
tab

##           M      Mean     Mean2 Prob(theta<0)
##  [1,] 1e+01 0.5376929 0.5508078     0.2065551
##  [2,] 2e+01 0.5319219 0.5572364     0.2314653
##  [3,] 3e+01 0.5319216 0.5572353     0.2345890
##  [4,] 4e+01 0.5319214 0.5572344     0.2355990
##  [5,] 5e+01 0.5319212 0.5572338     0.2360507
##  [6,] 6e+01 0.5319211 0.5572334     0.2362914
##  [7,] 7e+01 0.5319211 0.5572331     0.2364347
##  [8,] 8e+01 0.5319210 0.5572328     0.2365270
##  [9,] 9e+01 0.5319210 0.5572326     0.2365898
## [10,] 1e+02 0.5319209 0.5572325     0.2366346
## [11,] 1e+03 0.5319206 0.5572311     0.2368208
## [12,] 1e+04 0.5319206 0.5572309     0.2368227
## [13,] 1e+05 0.5319206 0.5572309     0.2368227
## [14,] 1e+06 0.5319206 0.5572309     0.2368227

PROBLEM 2: Monte Carlo integration

In this example, we want to compute \(Pr(\theta>0.7)\) when \(\theta \sim N(0,1)\) by using Monte Carlo integration based on \(M\) draws from \(N(0,1)\). More precisely, the integral we want to compute is \[ I = \int_\Re I(\theta>0.7)p(\theta)d\theta. \] The Monte Carlo integration approximation is \[ \widehat{I}_{MC} = \frac{1}{M} \sum_{i=1}^M I(\theta^{(i)}>0.7), \] where \(\theta^{(1)},\ldots,\theta^{(M)}\) is a random sample of \(N(0,1)\). It is easy to see, based on frequentist asymptotic results, that \[ \widehat{I}_{MC} \longrightarrow I \qquad \mbox{when} \ \ \ M \longrightarrow \infty. \] An important feature of Monte Carlo methods is the computation of the Monte Carlo error \[ V(\widehat{I}_{MC}) \approx \frac{1}{M^2} \sum_{i=1}^M (\theta^{(i)}-\widehat{I}_{MC})^2. \] The figure below illustrates the above asymptotic result by showing the improvement of the MC approximation as \(M\) increases from \(500\) to \(10,000\). The MC error descreases as the sample increases.

set.seed(314159)
M = 10000
theta.draw = rnorm(M,0,1)
I.mc       = cumsum(ifelse(theta.draw>0.7,1,0))/1:M
mc.error   = sqrt(mean(theta.draw-I.mc)^2/(1:M))
bands = cbind(I.mc-2*mc.error,I.mc,I.mc+2*mc.error)

M0 = 500
par(mfrow=c(1,1))
plot(M0:M,bands[M0:M,2],type="l",xlab="Monte Carlo size",
     ylab="MC approximation",ylim=c(0.2,0.3),axes=FALSE)
axis(2);axis(1,at=seq(M0,M,by=1000))
lines(M0:M,bands[M0:M,1],col=2)
lines(M0:M,bands[M0:M,3],col=2)
abline(h=1-pnorm(0.7),lty=2,col=4)

Replication exercise

In this replication exercise, commonly known as a Monte Carlo exercise, we see the decreasing size of the Monte Carlo error.

set.seed(314159)
par(mfrow=c(1,1))
plot(M0:M,I.mc[M0:M],type="l",xlab="Monte Carlo size",
     ylab="MC approximation",ylim=c(0.2,0.3),col=grey(0.8),
     axes=FALSE)
axis(2);axis(1,at=seq(M0,M,by=1000))
for (i in 1:100){
  theta.draw = rnorm(M,0,1)
  I.mc = cumsum(ifelse(theta.draw>0.7,1,0))/(1:M)
  lines(M0:M,I.mc[M0:M],col=grey(0.8))
}
abline(h=1-pnorm(0.7),lwd=3)
title("100 replications")

PROBLEM 3: MCIS

Monte Carlo integration (or simulation) via importance sampling is a very popular strategy when one is not able to directly sample from the target distribution, \(p(\theta)\). The importance sampling samples from a proposal distribution, \(q(\theta)\).

Here, we want to compute \(Pr(\theta>0.7)\) when \(\theta \sim N(0,1)\) using Monte Carlo integration based on \(M\) draws from the proposal distribution \(U(-20,20)\).

M = 10000
theta.draw = runif(M,-20,20)
I.mcis = mean(ifelse(theta.draw>0.7,1,0)*dnorm(theta.draw)/(1/40))
c(I.mc[M],I.mcis)

## [1] 0.2421000 0.2469929

PROBLEM 1 (CONTINUING)

Computing \(Pr(\theta<0)\) by Monte Carlo via Importance Sampling (MCIS) where the proposal density is \(N(0,4)\).

M = 100000
theta.draw  = rnorm(M,0,2)
denominator = exp(-0.5*(theta.draw-1)^2)*(1+0.25*theta.draw^2)^(-2.5)/dnorm(theta.draw,0,2)
numerator   = ifelse(theta.draw<0,1,0)*denominator
I.mc = mean(numerator)/mean(denominator)

paste("Exact value = ",round(Prob,5),sep="")

## [1] "Exact value = 0.23682"

paste("MCIS approximation = ",round(I.mc,5),sep="")

## [1] "MCIS approximation = 0.23632"

PROBLEM 4: Rejection method

Approximating the constant \(\pi\) via rejection method.

set.seed(3141593)
M = 10000
x = runif(M)
y = runif(M)
ind = order(x)
x=x[ind]
y=y[ind]
cond = (x^2+y^2)<1
plot(x,y,pch=16,cex=0.5)
points(x[cond],y[cond],pch=16,cex=0.5,col=2)
lines(x,sqrt(1-x^2))

pi.mc = 4*sum(cond)/M
pi.mc

## [1] 3.1312

PROBLEM 5: SIR

Sampling from the target distribution, \(N(0,1)\), based on draws from the proposal density, \(Cauchy(0,1)\) via Sampling Importance Resampling (SIR) and approximating the following few summaries:

• \(Pr(\theta<1)\)

• \(E(|\theta|)\)

• 2.5% quantile

Recall that the SIR algorithm works as follows:

1. Sample (large) \(M\) draws from a proposal density \(q(\theta)\): \(\{\widetilde{\theta}^{(1)},\ldots,\widetilde{\theta}^{(M)}\}\);

2. Evalute the weights \(\omega^{(i)} \propto \frac{p(\widetilde{\theta}^{(i)})}{q(\widetilde{\theta}^{(i)})}\);

3. Resample \(N \leq M\) from the discrete set \(\{\widetilde{\theta}^{(1)},\ldots,\widetilde{\theta}^{(M)}\}\) with probabilities \(\{\omega^{(1)},\ldots,\omega^{(M)}\}\).

4. The final sampled set, say \(\{\theta^{(1)},\ldots,\theta^{(N)}\}\) approximates a sample from \(p(\theta)\).

More details about Monte Carlo integration, MC simulation, SIR and other Monte Carlo schemes, can be found in my course notes here.

M = 100000
N = 50000
theta = rt(M,df=1)
w = exp(-theta^2/2)*(1+theta^2)
theta.new = sample(theta,size=N,replace=TRUE,prob=w)

par(mfrow=c(1,2))
hist(theta[abs(theta)<15],breaks=seq(-15,15,length=40),
prob=TRUE, xlab=expression(theta),main="")
title("Draws from proposal: Cauchy(0,1)")
legend("topleft",legend=paste("Tails = ",
100*round(1-sum(abs(theta)<15)/M,4),"%",sep=""),bty="n")

hist(theta.new,prob=TRUE,xlab=expression(theta),main="")
title("Draws from target: N(0,1)")
thetas = seq(-5,5,length=100)
lines(thetas,dnorm(thetas),col=2,lwd=3)

tab = cbind(c(pnorm(1),mean(theta.new<1)),
c(sqrt(2/pi),mean(abs(theta.new))),
c(qnorm(0.025),quantile(theta.new,0.025)))
rownames(tab) = c("Exact","SIR")
colnames(tab) = 
c("Pr(theta<1)","E(|theta|)","2.5% quant")
round(tab,4)

##       Pr(theta<1) E(|theta|) 2.5% quant
## Exact      0.8413     0.7979    -1.9600
## SIR        0.8415     0.7984    -1.9525