5 Week 10: Random walk Metropolis and applications of the MH algorithm

Let $X$ have pdf $f(x)$, where $f(x)=0$ for . Then multiplicative RWM essentially transforms the problem onto the log-scale as follows. Propose

$$\log Y\sim N(\log x,{\sigma }^2)\mathit{\hspace{1em}\hspace{1em}}\text{as }n\to \mathrm{\infty }.$$

Or equivalently, if $Z=\log X$ and $z=\log x$, propose

$$W\sim N(z,{\sigma }^2)\mathit{\hspace{1em}\hspace{1em}}\text{as }n\to \mathrm{\infty }.$$

Let $f_Z(\cdot )$ denote the pdf of $Z$, then

	$f_Z(z)$	$=$	$\left|{\displaystyle \frac{dx}{dz}}\right|f(x=e^z)$
		$=$	$e^z\pi (e^z)$

since $x=e^z$. Therefore for multiplicative RWM, the acceptance probability is

$$\min \{1,\frac{f_Z(W)}{f_Z(z)}\}=\min \{1,\frac{f_Z(\log Y)}{f_Z(\log x)}\}=\min \{1,\frac{Yf(Y)}{xf(x)}\}.$$

(The first term is the result of doing RWM upon the log scale, with the other two terms being the result of the transformation.)

We demonstrate the benefits of multiplicative RWM with

For both standard RWM and multiplicative RWM I started the MCMC algorithm in the tail of $X$ with $X_0=100$. (Note that $P(X>100)=1/1030301$.) Trace plots of the first 1000 iterations of each algorithm are presented below. We can see that the multiplicative RWM quickly converges to the centre of probability mass around 0, whereas the RWM algorithm shows random walk behaviour around 100. The RWM algorithm does eventually get close to 0 but it takes quite a few iterations. The RWM algorithm is correct but performs poorly for heavy tailed densities (probability densities where $f(x)e^{\theta x}\to \mathrm{\infty }$ as $x\to \mathrm{\infty }$ for all $\theta >0$). This is, because the RWM can get stuck in tails of the distribution and wanders apparently aimlessly for a long time before returning to the area of high probability mass. Since for large $x$, if $Y=x+ϵ$, where $ϵ$ is small compared to $x$, $f(Y)/f(x)\approx 1$. Thus most moves are accepted, hence the random walk behaviour.

Unnumbered Figure: Link

RWM with $\sigma =2$

Unnumbered Figure: Link

Multiplicative RWM with $\sigma =1$.

In addition, both algorithms were run started at $X_0=1$ for $100000$ iterations. Both algorithms provide good estimation of $f(\cdot )$ in the range $[0,10]$. However, if we consider large values of $x$, we can see that multiplicative RWM is better. In 100000 iterations, on average there will be 75 observations greater than 10 and 11 observations greater than 20. For the RWM algorithm there was only 23 observations greater than 10 and none greater than 20 (the maximum value was 13.5). On the other hand for multiplicative RWM there were 67 observations greater than 10 and 17 observations greater than 20. The maximum value for multiplicative RWM was 41 and under $f$, the probability of getting a value greater than 41 is $1/74088$. Thus the RWM algorithm is less inclined to enter the tail of $f$ than the multiplicative RWM algorithm. However, as we saw above when the RWM algorithm does enter the tail of $\pi$ it spends a long time in the tail before returning to the area of high probability mass.

The standard RWM algorithm proposes a new value centered on the current value. However, it makes sense, where possible, to make use of knowledge of $f$. In particular, we want to propose more values $Y$ where $f$ is large than where $f$ is small. Therefore if we know $f(x)$, and in particular, its derivative we can use this to concentrate the proposal in areas of the distribution with high probability mass. The MALA algorithm (Metropolis adjusted Langevin algorithm) is as follows.

1. 1.

   Propose a new value $y$, by setting

   $$y=x+\frac{{\sigma }^2}{2}\frac{d}{dx}\log f(x)+N(0,{\sigma }^2),$$

   where $X_t=x$.

2. 2.

   Accept the proposed value with probability

   $$\min \{1,\frac{f(y)\exp \left(-\frac{1}{2{\sigma }^2}{\left\{x-y-\frac{{\sigma }^2}{2}\frac{d}{dx}\log f(Y)\right\}}^2\right)}{f(x)\exp \left(-\frac{1}{2{\sigma }^2}{\left\{y-x-\frac{{\sigma }^2}{2}\frac{d}{dx}\log f(x)\right\}}^2\right)}\}.$$

3. 3.

   If the proposed value is accepted set $X_{t+1}=Y$, otherwise set $X_{t+1}=x$.

The key point to note is that the proposal is not symmetric, ie. $q(x,y)\ne q(y,x)$, and therefore this needs to be taken account of in the acceptance probability. The class of models to which MALA can be applied is smaller than RWM and using MALA is more involved as we need the derivatives of $f(\cdot )$. However, where MALA can be used it is preferable to RWM. For high-dimensional problems MALA is seen to have some very nice properties, see [4]. In particular, if $d$, the number of dimensions, is large, the optimal choice of ${\sigma }^2$ for MALA is $O(d^{-1/3})$ compared with $O(d^{-1})$ for RWM. Thus potentially offering substantial improvements.

Suppose that for $i=1,2,\mathrm{\dots },n$, $Y_i\sim \mathrm{Po}(\lambda )$. The likelihood is

$$\pi (𝐲|\lambda )\propto {\lambda }^{{\sum }_{i=1}^n{y}_i}{e}^{-n\lambda }.$$

The conjugate prior for $\lambda$ is a Gamma($\alpha ,\beta$) distribution $\pi (\lambda )\propto {\lambda }^{\alpha -1}{e}^{-\beta \lambda }$.

But gamma priors are often extremely unsatisfactory for expressing ignorance. (e.g. just plot the graph with $\alpha =1/2$.)

Instead use a half-Cauchy prior ${\pi }_0(\lambda )\propto 1/(1+{\beta }^2{\lambda }^2)(\lambda \ge 0)$, for some $\beta$ (e.g. $\beta =1$). The posterior is

$$\pi (\lambda |𝐲)\propto \frac{{\lambda }^{{\sum }_{i=1}^n{y}_i}{e}^{-n\lambda }}{1+{\beta }^2{\lambda }^2}.$$

Perhaps the easiest MH algorithm in this case would be an independence sampler with ${\lambda }^{*}\sim \text{Gam}(1+{\sum }_{i=1}^n{y}_i,n)$, i.e.

$$q(\lambda )\propto {\lambda }^{{\sum }_{i=1}^n{y}_i}{e}^{-n\lambda }.$$

Since in this case

$$\frac{\pi (\lambda |𝐲)}{q(\lambda )}\propto \frac{1}{1+{\beta }^2{\lambda }^2}\le 1.$$

Alternatively a random walk Metropolis algorithm with

$${\lambda }^{*}\sim N({\lambda }_{curr},{2.4}^2\times \frac{1}{n^2}\left(1+\sum _{i=1}^n{y}_i\right))$$

should also work well.
Why the above choice of proposal variance?
The variance of $\text{Gam}(1+\sum y_i,n)$ is $\frac{1}{n^2}(1+\sum y_i)$. (This is the posterior with the improper, flat prior $\pi (\lambda )\propto 1$ $(\lambda >0)$.) Scale by a factor of ${2.4}^2$ as we know this is the optimal proposal for a $N(0,1)$. Note that if $n$ is large by the central limit theorem $\text{Gam}(1+\sum y_i,n)\approx N(\{1+\sum y_i\}/n,\{1+\sum y_i\}/n^2)$.

Here are the traceplots for the suggested algorithms applied to the following data $𝐲=(6,8,8,3,6,3,6,4)$ (so $n=8,$ and $\sum y_i=44$), and starting with $\lambda =5.5$. (The MLE of a Poisson is $\widehat{\lambda }=\sum y_i/n=44/8=5.5$.)

Unnumbered Figure: Link

Independence Sampler

Unnumbered Figure: Link

RWM

Mixing seems reasonable. Acceptance rates are $0.864$ and $0.425$ respectively and effective sample sizes are $611.5$ and $233.3$.

Estimates of the posterior median and $95\%$ credibility intervals for $\lambda$ are

Independence Sampler: 5.35 (3.83,7.11).

RWM: 5.36 (3.88,7.13).

Note:. In ”real life” you would not use MCMC on the above target; there is a simpler and more appropriate technique, rejection sampling, which you will have met in the Bayesian Statistics course.

In this Section we consider using the Metropolis Hastings algorithm for logistic regression models. We use the Caesarian section example previously studied in Section 3.5 (using a Probit regression and data augmentation Gibbs sampler).

Response: Occurrence or non-occurrence of infection following birth by Caesarian section.

Covariates:

$\bullet$ $x_1=1$ if Caesarian section was not planned and $x_1=0$ otherwise;

$\bullet$ $x_2=1$ if there is presence of one or more risk factors, such as diabetes or excessive weight, and $x_2=0$ otherwise;

$\bullet$ $x_3=1$ if antibiotics were given as prophylaxis and $x_3=0$ otherwise.

In total: $251$ births.

We have aggregated all binary responses in each covariate category to binomial responses. Thus,

$$Y_i\sim \mathrm{Binomial}(n_i,p_i)$$

where $i$ is an index for each covariate combination, $y_i$ is the number of infections and $n_i$ the total number of observations for the $i^{th}$ covariate combination.

Suppose that we assume the binary logistic regression model,

$$\text{logit}\left(p_i\right)=\log \left(\frac{p_i}{1-p_i}\right)={\eta }_i={𝐱}_i^T𝜷,$$

where $p_i$ denotes the probability of an infection for each individual in the $i^{th}$ covariate combination and ${𝐱}_i^T=(1,x_{1i},x_{2i},x_{3i})$.

	$f(𝐲|𝜷)$	$=$	${\displaystyle \prod _{i=1}^n}f(y_i|𝜷)\propto {\displaystyle \prod _{i=1}^n}{p}_i^{y_i}{(1-p_i)}^{n_i-y_i}$
		$=$	${\displaystyle \prod _{i=1}^n}{\left[{\displaystyle \frac{\exp (𝐱_i{}^T𝜷)}{1+\exp (𝐱_i{}^T𝜷)}}\right]}^{y_i}{\left[{\displaystyle \frac{1}{1+\exp (𝐱_i{}^T𝜷)}}\right]}^{n_i-y_i}$
		$=$	${\displaystyle \prod _{i=1}^n}{\displaystyle \frac{\exp (𝐱_i{}^T𝜷\cdot y_i)}{{[1+\exp (𝐱_i{}^T𝜷)]}^{n_i}}}.$

A Maximum likelihood analysis in R gives the following results.

$$ Estimate Std. Error z value Pr(>|z|)
(Intercept) -1.8926 0.4124 -4.590 4.44e-06 ***
noplan $$ 1.0720 0.4253 2.520 0.0117 *
factor $$ 2.0299 0.4552 4.459 8.23e-06 ***
antib $$ -3.2544 0.4813 -6.762 1.37e-11 ***
---
Sig. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Null deviance: 83.491 on 6 degrees of freedom
Residual deviance: 10.997 on 3 degrees of freedom
AIC: 36.178

Number of Fisher Scoring iterations: 4

The 10.997 value for the residual deviance is not consistent with a ${\chi }_3^2$ distribution, since we can compute $P({\chi }_3^2>10.997)$ as

> 1-pchisq(10.997,3) [1] 0.01174211

So there is a fairly strong warning sign that the model may not fit adequately. To keep things simple we will ignore this (could consider interactions or probit regression, as seen earlier).

A Bayesian analysis places a prior on $𝜷={({\beta }_0,{\beta }_1,{\beta }_2,{\beta }_3)}^T$, usually a multivariate normal prior $𝜷\sim N_4({𝝁}_0,{𝐂}_0)$. The posterior distribution for $𝜷$ is

	$\pi (𝜷|𝐲)$	$\propto$	$f(𝐲|𝜷){\pi }_0(𝜷)$
		$\propto$	$\mathrm{\hspace{0.17em}}\left\{{\displaystyle \prod _{i=1}^n}{\displaystyle \frac{\exp (𝐱_i{}^T𝜷\cdot y_i)}{{(1+\exp (𝐱_i{}^T𝜷))}^{n_i}}}\right\}\exp \left(-{\displaystyle \frac{{(𝜷-{𝝁}_0)}^T{𝐂}_0^{-1}(𝜷-{𝝁}_0)}{2}}\right).$

The posterior density is a complicated, non-linear function of $𝜷$. How can we construct an MCMC algorithm to sample from this distribution?

1. 1.

   Update $\beta$ jointly: Independence sampler with multivariate normal (or better, $t_5$) candidate generator with mean equal to posterior mode and covariance matrix equal to inverse curvature at mode.
   $\bullet$ Newton-Raphson for $\pi (𝜷|𝐲)$;
   $\bullet$ relies on asymptotic normality of posterior distribution;
   $\bullet$ Multivariate normal may rarely propose values in the tails of the posterior, hence the preference for a $t$-distribution.

2. 2.

   Update $\beta$ jointly: Multivariate Gaussian random walk candidate generator, with covariance matrix equal to the inverse curvature at posterior mode times a (scalar) tuning factor, or some other pre-estimated covariance matrix
   $\bullet$ Needs tuning and pre-calculation of covariance estimate.

3. 3.

   Update $\beta$ component-wise via MH: Simple random walk proposal updating one $\beta$ component at a time.
   $\bullet$ Needs tuning;
   $\bullet$ problems with collinearity.

The code to implement algorithms 1, 2 and 3 are stored in files logisticIP.R, logisticMRWM.R and logisticRWM.R. The code to run the analysis is given in the file logisticRUN.R.
The output from running these algorithms for 5000 iterations is plotted in the files “logistic1.pdf”, “logistic2.pdf” and “logistic3.pdf”.

Estimates (posterior mean, standard deviation and posterior probability of exceeding 0) using the first algorithm (independence sampler – 5,000 samples after burn-in of 500) are given in the following table:

Coefficients:	mean	Std	probability
(Intercept)	-1.9544	0.4228	0
noplan	1.1071	0.4229	0.9968
factor	2.0955	0.467	1
antib	-3.3322	0.4867	0

The overall acceptance rate of this run was 87.6% which is very good for the independence sampler. The multivariate random walk algorithm works well when the covariance matrix of the proposal distribution reflects correlations between the parameters. For both the independence sampler and the multivariate random walk algorithms the same mechanism is used to find an appropriate covariance matrix for the proposal distribution. The componentwise (univariate) random walk Metropolis performs more poorly with slowly decaying correlations in the acf plot. However, the componentwise random walk Metropolis is run without tuning of the covariance matrix which is important for successful implementation of the other algorithms.

This example is taken from [1] and gives a non-linear model for oxygen demand in a biochemical experiment against incubation time.

Model:

$$y={\beta }_1(1-\exp (-{\beta }_{2}x))+ϵ$$

where $y$ is biochemical oxygen demand in $mg/l$, $x$ is incubation time in days and $ϵ\sim N(0,{\sigma }^2)$.

The data with predictive curve $y=213.0894(1-\exp (-0.5472x))$ is given below. This is the least squares estimated curve.

Unnumbered Figure: Link

Data source: www.itl.nist.gov/div898/strd/nls/data/boxbod.shtml

Random walk Metropolis was used to obtain a sample from the posterior distribution of $𝜽=({\beta }_1,{\beta }_2,\sigma )$ with a multivariate Gaussian proposal. Tuning of the proposal variance found that $1.39V$ worked well, where $1.39=2.4/\sqrt{3}$ and

is computed on the basis of pilot runs.

The estimated posterior means (standard deviations) of the parameters are $207(25.3),0.886(0.818)$ and $30.0(17.2)$ for ${\beta }_1$, ${\beta }_2$ and $\sigma$, respectively, based on a sample of size 10000. (A burn-in of 1000 iterations was used.) Finally, in the figure below ${\beta }_2$ is plotted against ${\beta }_1$, which clearly shows a non-linear relationship between the regression parameters.

Unnumbered Figure: Link

The first few lines in lab5code.r read in the data and create a response vector y and the design matrix X. Make sure you understand what the code is doing. Look at the data. Why do we subtract 1 when creating the second column of X?

There are three covariates: habitat, latitude, and elevation, and so there are four components to our parameter vector $𝜷$ - why?

We will use independent Gaussian priors for each of the four $\beta$s, each centered on $0$ and with the same variance, the scalar value V.prior. The function log.prior() should take the covariate vector beta and the scalar variance V.prior and should return the log of the prior density for beta. Work out the formula for the log prior (up to an arbitrary constant independent of beta) and so fill in the template for this function.

We wish to employ vague but proper priors so V.prior should be large but finite: a value of 100 should suffice in this instance.

We will be using the following model

$$Y_i\sim Po({\lambda }_i),{\lambda }_i=\exp ({𝐱}_i^t𝜷).$$

Recall that the density of a single Poisson random variable $Y_i$ which has mean ${\lambda }_i$ is

$$f(y_i)=\frac{{\lambda }_i^{y_i}{e}^{-{\lambda }_i}}{y_i!}.$$

Use vector arithmetic to fill in the template for the function poisson.log.like() which returns the log-likelihood of the parameter vector beta given data vector y and covariate matrix X. Hint: your very first step should be loglam <- X %*% beta. This gives a vector of the $\log (\lambda )$ values.

The inputs and outputs of the function RWM() are detailed in the comments above it. Fill in the blanks as instructed by the comments in the function.

Note V.prop is the variance of the (multivariate) Gaussian jump proposal in the RWM for beta. Since beta is a 4-component vector, V.prop must be a 4x4 matrix.

Initially you will tune the algorithm by choosing a scaling and a variance matrix. All tuning runs should be for 10000 iterations and should start at $\beta =(0,0,0,0)$. Set your plot window to $2\times 2$ with the command par(mfrow=c(2,2)).

1. 1.

   Perform a run with V.prop=diag(rep(1,4)) (why?) and lambda.prop=0.001. Create an MCMC object from this run and look at the trace plots. What do you notice?

2. 2.

   Find the variance matrix from this run using the command V=var(B) where B is the output from the previous run. Change your scale parameter, lambda.prop, to 1 (this is approximately $2.38/\sqrt{4}$) and try a new run with V.prop set to V. What do the traceplots show?

3. 3.

   Repeat 2. until all four components are mixing well. What is your final variance matrix? What is the correlation matrix? Use these to explain the behaviour of the first run.

4. 4.

   Look at the effective sample size for this run (removing burn in by slicing). Can you increase it by changing lambda.prop?

5. 5.

   Use slicing on your final run to chop the first few rows from your output matrix (burn-in). Find the posterior median estimate and a $95\%$ CI for each parameter and draw kernel density estimates or histograms from the MCMC output.