Inferring the error variance in Metropolis-Hastings MCMC

One of the great joys of working with two talented post-docs in the research group – Mike Stout and Mudassar Iqbal – as well as a great collaboration with Theodore Kypraios, is that they are often one step ahead of me and I am playing catch-up. Recently, Theo has discussed with them how to estimate the error variance associated with the data used in Metropolis-Hastings MCMC simulations.

The starting point, usually, is that we have some data, let us say y_i for i=1, \cdots, n, and a model – usually, in our case, a dynamical system – which we are trying to fit to the data. For any given set of parameters \theta, our model will provide estimates for the data points that we will call \hat y_i. Now, assuming uniform Gaussian errors, our likelihood function L(\theta) looks like:

L(\theta) = \prod_{i=1}^n \frac{1}{\sqrt{2 \pi\sigma^2}}e^{-\frac{1}{2}(\frac{y_i - \hat y_i}{\sigma})^2}

where \sigma^2 is the error variance associated with the data. Now, when I first started using MCMC, I naively thought that we could use values for \sigma^2 provided by our experimental collaborators, and so we could use different values of \sigma^2 according to how confident our collaborators were in the measurements, equipment etc. What I found in practice was that these values rarely worked (in terms of convergence of the Markov chain) and we have had to make up error variances using trial and error.

So I was delighted when I heard that Theo had briefed both Mike and Mudassar about a method for estimating the error variance as part of the MCMC. Since I have not tried it before, I thought I would give it a go. I am posting the theory and some of my simulations, which are helpful results.

Theory

The theory behind estimating \sigma^2 is as follows. First, set

\tau = \frac{1}{\sigma^2}

We can then re-write the likelihood, now for the model parameters \theta and also the unknown value \tau, as

L(\bf{\theta}, \tau) = \frac{\tau^{(n/2)}}{\sqrt{2 \pi}^n}e^{-\frac{\tau}{2}\sum_{i=1}^n(y_i - \hat y_i)^2}

Now observe that this has the functional form of a Gamma distribution for \tau, as the p.d.f. for a Gamma distribution is given by:

f(x; \alpha, \beta) = \frac{\beta^\alpha}{\Gamma(\alpha)}x^{\alpha-1}e^{-\beta x}

So if we set a prior distribution for \tau as a Gamma distribution with parameters \alpha and \beta, then the conditional posterior distribution for \tau is given by:

p(\tau | \theta) \propto \tau^{(n/2)+ \alpha - 1}e^{-\tau(\frac{1}{2}\sum_{i=1}^n(y_i - \hat y_i)^2+\beta)}

We observe that this is itself a Gamma distribution, with parameters \alpha \prime = \alpha + n/2 and \beta \prime = \beta + \frac{1}{2} \sum_{i=1}^n (y_i - \hat y_i)^2. Thus the parameter \tau can be sampled with a Gibbs step as part of the MCMC simulation (usually using Metropolis-Hastings steps for the other parameters).

Simulations

The simulations I have run are with a toy model that I use a great deal for teaching. Consider a constitutively-expressed protein that is produced at constant rate k and degrades (or dilutes) at constant rate \gamma per protein. A differential equation for protein concentration P is given by:

\frac{dP}{dt} = k - \gamma P

This ODE has the closed form solution:

P = \frac{k}{\gamma} + (P_0 - \frac{k}{\gamma}) e^{-\gamma t}

where P_0 is the concentration of protein at t=0. For the purposes of MCMC estimation, mixing is improved by setting P_1 = \frac{k}{\gamma} so that the closed form solution is:

P = P_1 + (P_0 - P_1) e^{-\gamma t}

Some data I have used for teaching purposes comes from the paper Kim, J.M. et al. 2006. Thermal injury induces heat shock protein in the optic nerve head in vivo. Investigative ophthalmology and visual science 47: 4888-94. The data is quantitative Western blots of Hsp70 in the optic nerve of rats, as induced by laser damage. (Apologies for the unpleasantness of the experiment):

Time / hours Protein / au
3 1100
6 1400
12 1700
18 2100
24 2150

The aim is to use a Metropolis-Hastings MCMC, together with a Gibbs step for the \tau parameter, to fit the data. The issue that immediately arises is how to set the parameters \alpha and \beta. This may seem arbitrary, but it is already better than choosing a value for \sigma^2, as the Gamma distribution will exploring of that parameter. For my first simulation, I thought that \sigma = 100 would be sensible (this turned out to be a remarkably good choice, as we will see). So I set \alpha = 0.01 and \beta = 100 and lo and behold, the whole MCMC worked beautifully. (Incidentally, I used independent Gaussian proposals for the other three parameters, with standard deviations of 100 for the P_0 and P_1 and standard deviation of 0.01 for \gamma. These parameters were forced to be positive – Darren Wilkinson has an excellent post on doing that correctly. Use of log-normal proposals in this case leads to very poor mixing, with the chain taking some large excursions for the P_1 and \gamma parameters).

histograms

The median parameter values are P_0 = 786, P_1 = 2526, \gamma = 0.0686 and \tau = 0.000122. The latter corresponds to \sigma = 90.6. With these values, we can see a good fit to the data: below are plotted the data points (in red), the best fit (with median parameter values) in blue, and model predictions from a random sample of 50 parameter sets from the posterior distribution in black.

fit

Considerations

However, some questions obviously arise: how sensitive is this procedure to choices of \alpha and \beta? I will confess: I use Bayesian approaches fairly reluctantly, being more comfortable with classical frequentist statistics. What I like about Bayesian approaches are firstly the description of unknown parameters with a probability distribution, and secondly the availability of highly effective computer algorithms (i.e. MCMC). What makes me uncomfortable is the potential for introducing bias through the prior distributions. So I have carried out some investigations with different values of \alpha and \beta. In particular, I wanted to know: (i) what happens if I keep the mean (equal to \alpha / \beta) the same but vary the parameters? (ii) what happens if I vary the mean of the distribution? The table below summarizes positive results:

alpha beta P0 P1 gamma sigma
0.01 100 786 2526 0.0686 90.6
1 10000 747 2428 0.0795 98.0
0.0001 1 797 2533 0.0681 96.3
0.1 10 822 2603 0.0623 97.9
0.001 1000 760 2455 0.0758 94.8
1 1 792 2539 0.0676 64.9
0 0 805 2565 0.0653 98.3

As you can see (please ignore the last line for now), the results are robust to a very wide range of \alpha and \beta, even producing a good estimate for \sigma when that estimate is a long way from the mean of the prior distribution. But then we can make the following observation. Consider the sum of squares for a ‘best-fit’ model, for example using the parameters for the first row (this is 12748). So as long as \alpha \ll n/2 and \beta \ll 12748/2, the prior will introduce very little bias. But if you try to use values of \alpha and especially \beta very much larger than an estimated sum of squares from well-fitted model parameters, then things might go wrong. For example, when I set \alpha = 1 and \beta = 10^6 then my MCMC did not converge properly.

This leads to my final point, and the final row in the table. Would it be possible to remove prior bias altogether? If you look at the marginal posterior for \tau, we observe that if we set \alpha = \beta = 0, we obtain a Gamma distribution, whose mean is precisely the error variance, as, in this case,

\frac{\beta \prime}{\alpha \prime} = \frac{\sum_{i=1}^n(y_i - \hat y_i)^2}{n}

The algorithm should work perfectly well sampling from this Gamma distribution, and indeed it does, producing comparable results to when an informative prior is used.

Conclusions

In summary, I am happy to conclude that this method is good for estimating error variance. Clear advantages are:

  1. It is simple to implement and fairly fast to run – adding a Gibbs step is no big deal.
  2. It is clearly preferable to making up a fixed number for the error variance – which was what we were doing before.
  3. The prior parameters allow you to make use of information you might have from experimental collaborators on likely errors in the data.
  4. The level of bias from the priors is relatively low, and can be eliminated altogether.
Advertisements

One thought on “Inferring the error variance in Metropolis-Hastings MCMC

  1. Pingback: Heat Shock Proteins and the Brain: Implications for Neurodegenerative Diseases and Neuroprotection · My-infoWorld.NET

Leave a Reply

Fill in your details below or click an icon to log in:

WordPress.com Logo

You are commenting using your WordPress.com account. Log Out / Change )

Twitter picture

You are commenting using your Twitter account. Log Out / Change )

Facebook photo

You are commenting using your Facebook account. Log Out / Change )

Google+ photo

You are commenting using your Google+ account. Log Out / Change )

Connecting to %s