Random-walk Metropolis Algorithm

Nonlinear regression example

1 Unnormalized posterior density

The following example was taken from Chapter 6 (Metropolis-Hastings algorithms) of our book Gamerman and Lopes (2006) Markov Chain Monte Carlo: Stochastic Simulation for Bayesian Inference.

The unnormalized posterior (prior times likelihood) is given by \[ p(\theta|\mbox{data}) \propto p_N(\theta,0,10) \times \prod_{i=1}^n p_N(y_i,50+170x_i/(\theta+x_i),126), \] where \(\mbox{data}=(x,y)\), \(x=(0.02,0.02,0.06,0.06,0.11,0.11,0.22,0.22,0.56,0.56,1.10,1.10)\) and \(y=(76,47,97,107,123,139,159,152,191,201,207,200)\) are non-linearly related with mean function and residual variance given by, \[ E(y|x,\theta) = 50+\frac{170x}{\theta+x} \ \ \ \mbox{and} \ \ \ V(y|x,\theta)=126, \] respectively. The only unknown quantity/parameter is \(\theta\).

posterior = function(theta){
  dnorm(theta,0,10)*prod(dnorm(y,50+170*x/(theta+x),sqrt(126)))
}
logposterior = function(theta){
  dnorm(theta,0,10,log=TRUE)+sum(dnorm(y,50+170*x/(theta+x),sqrt(126),log=TRUE))
}
x = c(0.02,0.02,0.06,0.06,0.11,0.11,0.22,0.22,0.56,0.56,1.10,1.10)
y = c(76,47,97,107,123,139,159,152,191,201,207,200)

1.1 Pointwise evaluation

thetas = seq(0.08,0.20,length=1000)
post = rep(0,1000)
for (i in 1:1000)
  post[i] = posterior(thetas[i])

par(mfrow=c(1,1))
plot(thetas,post,xlab=expression(theta),ylab="Unnormalized posterior",type="l")

2 Posterior draws via random walk Metropolis

th    = 0.08
burnin = 1000
M = 10000
niter = burnin+M
th.d  = rep(0,niter)
for (iter in 1:niter){
  th.star = rnorm(1,th,0.01)
  logalpha = min(0,logposterior(th.star)-logposterior(th))
  if (log(runif(1))<logalpha){
    th = th.star
  }
  th.d[iter] = th
}
th.d = th.d[(burnin+1):niter]

2.1 Trace and ACF plots

par(mfrow=c(1,2))
ts.plot(th.d,xlab="Iterations",ylab=expression(theta))
acf(th.d,main="")

2.2 Histogram approximation to posterior

par(mfrow=c(1,2))
hist(th.d,prob=TRUE,main="",xlab=expression(theta),ylab="Posterior")
plot(thetas,post,xlab=expression(theta),ylab="Unnormalized posterior",type="l")

3 Computing \(p(y_{new}|x_{new},x_{1:n},y_{1:n})\)

Below we approximate via Monte Carlo the posterior preditives \[ p(y_{new}|x_{new},x_{1:n},y_{1:n})= \int_\Theta p(y_{new}|x_{new},\theta)p(\theta|x_{1:n},y_{1:n})d\theta, \] for a range of values of \(x_{new}\) that cover the observed ones in \(x_{1:n}\). In particular, the predictive mean function \[ E(Y|x_{new},x_{1:n},y_{1:n}) = \int_\Theta E(Y|x_{new},\theta)p(\theta|x_{1:n},y_{1:n})d\theta, \] can be approximated by \[ \frac{1}{M} \sum_{i=1}^M E(Y|x_{new},\theta^{(i)}), \] where \(\{\theta^{(1)},\ldots,\theta^{(M)}\}\) are the MC draws from \(p(\theta|x_{1:n},y_{1:n})\). This is the solid red line in the figure below. The dashed red lines are the 5th and 95th percentiles of \(p(y_{new}|x_{new},x_{1:n},y_{1:n})\). For comparison we added the MC approximation to \[ E(Y|x_{new},{\tilde \theta}), \] where \({\tilde \theta}=E(\theta|x_{1:n},y_{1:n})\). This is the solid blue line in the figure below. This is an example of the results that, in general, \[ E(g(X)) \neq g(E(X)). \]

xx = seq(min(x),max(x),by=0.01)
nxx = length(xx)
fx = matrix(0,M,nxx)
for (i in 1:M)
  fx[i,] = 50+170*xx/(th.d[i]+xx)
fx.mean = 50+170*xx/(mean(th.d[i])+xx)

par(mfrow=c(1,1))
plot(x,y)
for (i in 1:M)
  lines(xx,fx[i,],col=grey(0.95))
lines(xx,apply(fx,2,mean),col=2)
lines(xx,apply(fx,2,quantile,0.05),col=2,lty=2)
lines(xx,apply(fx,2,quantile,0.95),col=2,lty=2)
lines(xx,fx.mean,col=4)
points(x,y,pch=16)