Monte Carlo integration

Simple linear regression

1 Normal, Cauchy and Laplace priors

A parameter \(\beta\) is standard normal, Cauchy (standard Student’s \(t\) with one degree of freedom) or standard double-exponential (aka Laplace) if \[\begin{eqnarray*} p(\beta) &=& \frac{1}{\sqrt{2\pi}}\exp\{-0.5\beta^2\}, \ \ \mbox{or} \\ p(\beta) &=& \frac{1}{\pi(1+\beta^2)}, \ \ \mbox{or} \\ p(\beta) &=& \frac12 \exp\{-|\beta|\}, \end{eqnarray*}\] respectively, for \(\beta \in \Re\).

betas = seq(-5,5,length=1000)
plot(betas,dnorm(betas),type="l",ylim=c(0,0.5),xlab=expression(beta),ylab="Density",lwd=2)
lines(betas,dt(betas,1),col=2,lwd=2)
lines(betas,0.5*exp(-abs(betas)),col=3,lwd=2)
legend("topleft",legend=c("N(0,1)","t(1)","DE(1)"),col=1:3,lwd=2,lty=1,bty="n")

2 A simple linear regression example

set.seed(12345)
n    = 10
beta = 2
x    = rnorm(n)
y    = beta*x + rnorm(n)

plot(x,y)
abline(lm(y~x-1),col=2)
abline(h=0,lty=2)
abline(v=0,lty=2)

2.1 Maximum likelihood estimate

The maximum likelihood estimation of the simple Gaussian linear regression through the origin \[ y_i|x_i,\beta \sim N(\beta x_i,1), \] for \(i=1,\ldots,n\), is straightforward \[ {\widehat \beta}_{mle} = \frac{\sum_{i=1}^n y_ix_i}{\sum_{i=1}^n x_i^2}, \] such that \({\widehat \beta}_{mle}|\beta \sim N(\beta,v^2)\), where \(v^2=\frac{1}{\sum_{i=1}^n x_i^2}\).

beta.mle = sum(x*y)/sum(x^2)
se.beta  = 1/sqrt(sum(x^2))
c(beta.mle,se.beta)

## [1] 1.6650399 0.4037297

3 Bayesian inference via MCI

One of the main (computation) problems faced by applied Bayesian modelers/scientists is the computation of integrals, such as the normalizing constant \[ p(y|x) = \int_{-\infty}^\infty \left[\prod_{i=1}^n p_N(y_i|x_i\beta,1)\right]p(\beta)d\beta, \] or the posterior mean, \[ E(\beta|y,x) = \frac{\int_{-\infty}^\infty \beta \left[\prod_{i=1}^n p_N(y_i|x_i\beta,1)\right]p(\beta)d\beta}{p(y|x)}, \] or even \(Pr(\beta<0.75|y,x)\), \[ Pr(\beta<0.75|y,x) = \frac{\int_{-\infty}^\infty 1(\beta<0.75) \left[\prod_{i=1}^n p_N(y_i|x_i\beta,1)\right]p(\beta)d\beta}{p(y|x)}. \] Basically, Monte Carlo integration approximates the integral \[ I = \int h(\beta)g(\beta)d\beta = \int \left[\frac{h(\beta)g(\beta)}{q(\beta)}\right]q(\beta)d\beta, \] based on draws \(\{\beta^{(i)}\}_{i=1}^M\) from \(g(\beta)\) or based on draws \(\{{\tilde \beta}^{(i)}\}_{i=1}^M\) from \(q(\beta)\) (proposal/candidate/importance distribution), as \[\begin{eqnarray*} {\hat I}_1 &=& \frac{1}{M} \sum_{i=1}^M h(\beta^{(i)})\\ {\hat I}_2 &=& \frac{1}{M} \sum_{i=1}^M \frac{h({\tilde \beta}^{(i)})g({\tilde \beta}^{(i)})}{q({\tilde \beta}^{(i)})}. \end{eqnarray*}\] As \(M \longrightarrow \infty\), both \({\hat I}_1\) and \({\hat I}_2\) converge to \(I\).

Below we use a uniform distribution in the interval $(-5,5) as our proposal \(q(\beta)\).

M       = 10000
b.draws = runif(M,-5,5)
hist(b.draws,prob=TRUE,main="Draws from the proposal U(-5,5)",xlab=expression(beta))

3.1 Computing prior predictives

Here we compute \(p(y|x)\), the prior predictive density, based on each one of the three prior specifications of \(\beta\) (Gaussian, Student’s \(t\) and double-exponential).

like = rep(0,M)
for (i in 1:M)
  like[i] = prod(dnorm(y,x*b.draws[i],1))
k.n = mean(like*dnorm(b.draws)/(1/10))
k.t = mean(like*dt(b.draws,df=1)/(1/10))
k.l = mean(like*0.5*exp(-abs(b.draws))/(1/10))
c(k.n,k.t,k.l)

## [1] 4.475916e-07 3.650963e-07 4.087623e-07

The only one we can obtain in closed-form is under the normal prior, \(\beta \sim N(0,1)\): \[ p(y|x,\mbox{normal prior}) = p_N(y|0,xx'+I_n) \]

#install.packages("mvtnorm")
library("mvtnorm")
XtX = x%*%t(x)+diag(1,n)
k.n.exact = dmvnorm(y,rep(0,n),XtX,log=FALSE)
c(k.n.exact,k.n)

## [1] 4.484789e-07 4.475916e-07

3.2 Computing posterior means

E.n = mean(b.draws*like*dnorm(b.draws)/(1/10))/k.n
E.t = mean(b.draws*like*dt(b.draws,df=1)/(1/10))/k.t
E.l = mean(b.draws*like*0.5*exp(-abs(b.draws))/(1/10))/k.l
c(E.n,E.t,E.l)

## [1] 1.429081 1.516261 1.499867

The only one we can obtain in closed-form is under the normal prior, i.e. \[ E(\beta|y,x,\mbox{normal prior}) = \frac{\sum_{i=1}^n y_ix_i}{\sum_{i=1}^n x_i^2}. \]

E.n.exact = sum(y*x)/(sum(x^2)+1)
c(E.n.exact,E.n)

## [1] 1.431680 1.429081

3.3 Computing posterior tails

P.n = mean((b.draws<0.75)*like*dnorm(b.draws)/(1/10))/k.n
P.t = mean((b.draws<0.75)*like*dt(b.draws,df=1)/(1/10))/k.t
P.l = mean((b.draws<0.75)*like*0.5*exp(-abs(b.draws))/(1/10))/k.l
c(P.n,P.t,P.l)

## [1] 0.03483404 0.03008621 0.03181092

The only one we can obtain in closed-form is under the normal prior

P.n.exact = pnorm(0.75,E.n.exact,sqrt(1/(sum(x^2)+1)))
c(P.n.exact,P.n)

## [1] 0.03431329 0.03483404

3.4 Priors, likelihood and posteriors

likelihood = rep(0,1000)
for (i in 1:1000)
  likelihood[i] = prod(dnorm(y,x*betas[i],1))
post.n = dnorm(betas)*likelihood/k.n
post.t = dt(betas,df=1)*likelihood/k.t
post.l = 0.5*exp(-abs(betas))*likelihood/k.l

plot(betas,post.n,type="l",col=2,lwd=2,xlab=expression(beta),ylab="Density")
lines(betas,post.t,col=3,lwd=2)
lines(betas,post.l,col=4,lwd=2)
lines(betas,dnorm(betas),lty=2,col=2,lwd=2)
lines(betas,dt(betas,1),col=3,lty=2,lwd=2)
lines(betas,0.5*exp(-abs(betas)),col=4,lty=2,lwd=2)
lines(betas,dnorm(betas,beta.mle,se.beta),lwd=2)
legend("topleft",legend=c("N(0,1) prior","t(1) prior","DE(1) prior","Likelihood"),
       col=c(2:4,1),lwd=2,lty=c(2,2,2,1),bty="n")
abline(v=0.75,lty=2)       

3.5 Quality of MC approximation

M       = 100000
b.draws = runif(M,-5,5)
like = rep(0,M)
for (i in 1:M)
  like[i] = prod(dnorm(y,x*b.draws[i],1))
M1 = c(seq(100,1000,by=100),seq(2000,10000,by=1000),seq(20000,M,by=10000))
P.n = NULL
for (m in M1){
  k.n = mean(like[1:m]*dnorm(b.draws[1:m])/(1/10))
  P.n = c(P.n,mean((b.draws[1:m]<0.75)*like[1:m]*dnorm(b.draws[1:m])/(1/10))/k.n)
}
plot(M1/100,P.n,xlab="Monte Carlo size (x100)",ylab="Integral",type="b")
abline(h=P.n.exact,col=2)

4 Monte Carlo simulation via SIR

Draws from a candidate distribution \(q(\beta)\) can be reweighted and resampled in order to approximate draws from a target distribution \(p(\beta)\). This is the SIR scheme, one of the most popular MC sampling method for situations where dim(\(\beta\)) is small, say less than 10. SIR comprises three simple steps:

• Step 1: Draw \(\{{\tilde \beta}^{(i)}\}_{i=1}^M\) from the proposal/candidate/importance distribution \(q(\beta)\).
• Step 2: For \(i=1,\ldots,M\), compute (unnormalized) weights \[ \omega^{(i)} = \frac{p({\tilde \beta}^{(i)})}{q({\tilde \beta}^{(i)})}. \]
• Step 3: \(\{\beta^{(i)}\}_{i=1}^M\) from \(\{{\tilde \beta}^{(i)}\}_{i=1}^M\) with weights \(\{\omega^{(i)}\}_{i=1}^M\).

Notice that, one only need to be able to evaluate both \(p(\cdot)\) and \(q(\cdot)\) up to normalizing constants since when the weights are normalized these constants cancel out. Additionally, one needs to easily sample from \(q(\cdot)\).

4.1 SIR for linear regression

Usually, \(p(\beta)\) is proportional to the product of the prior \(\pi(\beta)\) and the likelihood \(L(\beta;\mbox{data})\), i.e. \[ p(\beta) = \kappa \pi(\beta)L(\beta;y,x), \] where \(k\) is the normalizing constant. If, in addition, one decides to use the prior as proposal (not a smart idea in most cases!), i.e. \(q(\beta) = \pi(\beta)\), it is easy to see that the SIR weights are \[ \omega^{(i)} = L({\tilde \beta}^{(i)};y,x). \] Let us now use SIR to sample from the posterior of \(\beta\) in our regression example when using the double exponential prior

4.1.1 Step 1 - Drawing from proposal

Here we will use the \(U(-5,5)\) as \(q(\beta)\).

M       = 20000
b.draws = runif(M,-5,5)

4.1.2 Step 2 - Computing weights

w = rep(0,M)
for (i in 1:M)
  w[i] = prod(dnorm(y,x*b.draws[i],1))*0.5*exp(-abs(b.draws[i]))/(1/10)

4.1.3 Step 3 - Resampling

beta.post = sample(b.draws,size=M,replace=TRUE,prob=w)
hist(beta.post,prob=TRUE,xlab=expression(beta),ylab="Density",main="")
lines(betas,post.l,col=2,lwd=2)

Now we can compare the MC integration with the MC simulation approximations to the three quantities of interest from the posterior \(p(\beta|y,x)\)

c(E.l,mean(beta.post))

## [1] 1.499867 1.501845

c(k.l,mean(w))

## [1] 4.087623e-07 4.133555e-07

c(P.l,mean(beta.post<0.75))

## [1] 0.03181092 0.03260000