2
$\begingroup$

I'm using Metropolis-Hastings MCMC to find the set of parameters $\theta$ of a model $M$ that best fits my experimental data $x$ (with noise), where it is possible to directly calculate $x$ (without noise) from $\theta$ but not to directly calculate $\theta$ from $x$. By sampling the (posterior) probability distribution of each $\theta$ being the best fit, which is $P(\theta|x)$, MCMC will also let me look at things like parameter correlations and sensitivity of each parameter to the fit. First of all, I'd like to confirm whether my basic knowledge about MCMC is correct.

The way the Metropolis-Hastings algorithm works is that if I know a function $f(\theta|x)$ that is proportional to the posterior, I can sample the posterior by starting at $\theta_0$ and moving with a proposal density $g(\theta_{i+1}|\theta_i)$ that suggests a candidate set of parameters given the current set of parameters. Each time a proposal is generated, I calculate the acceptance ratio $\alpha=f(\theta_{i+1}|x)/f(\theta_i|x)$ that determines whether to stay at $\theta_i$ or go to $\theta_{i+1}$. (I think $f$ is called the sampling density but I'm not sure.)

For the proposal density $g$, I'm using a Gaussian centered at $\theta_i$ with a fixed standard deviation $s$ that I specify, whish seems common from what I've read. The question is what to use for $f$. Since the goal is to minimize the error between $x$ and the calculated data $x_M(\theta)$, I have a least-squares fitness function, $F(\theta)=\frac{1}{N}\sum^N\left((x-x_M(\theta))^2/x^2\right)$, where $N$ is the number of points in $x$ and $x_M$. The reason it's $x^2$ in the denominator is to normalize the magnitude of $x$, which can be scaled differently depending on how one processes the data. Then I assume that the posterior is Gaussian, so for $f$ I use a Gaussian centered at $F=0$ with standard deviation $\sigma$, where $\sigma$ comes from the model. Therefore, the acceptance ratio is $\alpha=\mathrm{exp}\left(-\frac{1}{2\sigma^2}(F(\theta_{i+1})^2-F(\theta_i)^2)\right)$.

My question is, what is the meaning of $\sigma^2$ in this case? Does it describe the variance in the noise of the data (since $F$ depends partly on the variance in the data) or the variance in how close the parameters are to the best fit parameters (since $F$ also depends on this)? (Based on my least-squares equation, for $F=0$ there would have to be no noise in the data, and the fitted parameters would all have to be perfectly correct.) I'm thinking it's the latter case because it's basically a rewording of the posterior distribution, but in that case I can't wrap my head around how to model $\sigma^2$. I've been reading this page that towards the bottom talks about fitting $\sigma$ as a parameter, which confuses me even more.

$\endgroup$
  • $\begingroup$ $f$ should come from your model, you should not be tuning it for algorithmic performance of the MCMC. (If this is for sampling from a posterior distribution of a Bayesian model, which it seems to be as you talk about 'the posterior') $\endgroup$ – Juho Kokkala Aug 6 '16 at 4:39
  • $\begingroup$ Thanks for the comment. I wasn't tuning $f$ for performance but rather to get an accurate sampling from the posterior distribution. It sounds like this is not possible, and $f$ has to come from the model. I'll edit the question to ask about what $f$ can mean specifically for my case. $\endgroup$ – firescape Aug 8 '16 at 14:39
  • 1
    $\begingroup$ Changing $f$ changes the posterior distribution! $\endgroup$ – Juho Kokkala Aug 8 '16 at 15:12
  • $\begingroup$ How familiar are you with Bayesian inference? Note that MCMC (when used in the Bayesian context) is just a computational tool, while the issue here seems to be about specifying a model in the first place. $\endgroup$ – Juho Kokkala Aug 8 '16 at 16:17
  • $\begingroup$ After thinking about it some more, I need to incorporate the uncertainty of the observed data into the likelihood function, which is multiplied with the prior to get $f$. This is what has been done in the literature I've found. All the parameters of the likelihood function are a part of the model (including the $\sigma$ I was referring to in the case of a Gaussian likelihood). Anyway, you've answered the initial question, that $f$ should come from the model. $\endgroup$ – firescape Aug 8 '16 at 22:12
4
$\begingroup$

There are several papers that propose methods to adapt the Metropolis sampler in order to achieve certain optimality criterion. The most popular one consists of controlling the acceptance rate. Under certain regularity conditions of the target distribution, it has been shown that the Metropolis algorithm performs well if the acceptance rate is $\approx 0.44$ in dimension 1, and $\approx 0.234$ in higher dimensions:

Weak convergence and optimal scaling of random walk Metropolis algorithms (1997). A. Gelman, W. R. Gilks, and G. O. Roberts.

On the other hand, other self-adaptive algorithms that make use of running parallel chains have been developed in:

A general purpose sampling algorithm for continuous distributions (the t-walk) (2010). J. Andrés Christen and Colin Fox

Also implemented in the R package Rtwalk.

There is no general criterion to make efficient the Metropolis algorithm in every single scenario, but these two algorithms work well in non-pathological scenarios and low to moderate dimensions (<20).

$\endgroup$
  • 1
    $\begingroup$ Welcome to our site - thanks for taking the time to give full citations for your references, that's always appreciated. $\endgroup$ – Silverfish Aug 5 '16 at 23:02
  • $\begingroup$ This answers how to tune the variance of the proposal. But, as far as I understand the question, it seems that the standard deviation of the proposal $g(\cdot \mid \cdot)$ is fixed at $s$ and the OP is trying to select the standard deviation of something related to the target distribution $\endgroup$ – Juho Kokkala Aug 6 '16 at 4:44

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.