I'm coding a McMC algorithm for geophysical applications.

Using the Metropolis-Hastings scheme to accept/reject the proposed models is smth that i thought i completely understood, but i don't. To be clear: normally you have a current model $m$, you perturb one parameter (in my case velocity $v_i$) and in this way you get a trial model $m'$.

The acceptance probability for such a perturbation is:

$$\alpha(m'|m)=min\left[ 1, \frac{p(m')}{p(m)} \cdot \frac{p(d|m')}{p(d|m)} \cdot \frac{q(m|m')}{q(m|m)} \cdot |J|\right]$$

Where $p(m)$ are the priors, $p(d|m)$ are the likelihoods, $q(m|m')$ are the proposals and $J$ the jacobian.

So assuming that the priors are symmetric, and the jacobian is 1, I remain with:

$$\alpha(m'|m)=min\left[ 1,\frac{p(d|m')}{p(d|m)} \cdot \frac{q(m|m')}{q(m|m)} \right]$$

The likelihood ratio is not a problem, so now i wanna deal only with the proposal ratio:

Supposing that the way i make a perturbation is choosing one parameter, and perturb it with gaussian probability, so that the new velocity value will be $$v'_i=v_i + u \cdot \sigma$$ ($u$ is normal distributed with mean=0 and var.=1)

This perturbation can be described by a proposal in the form:

$$q(v'|v)= \frac{1}{\sigma \sqrt{2 \pi}} \cdot \exp \left[- \frac{(v'-v)^2}{2 \sigma ^2} \right]$$

It's easy to see that for this kind of perturbation the proposals are symmetric so the acceptance probability is just the ratio of the likelihoods.

Now the questions:

1) If instead of perturbing one parameter I want to change 3 of them at the same time, what happens to my proposals, and so to my acceptance probability $\alpha$? How should I compute them? ($v_1, v_2, v_3$ are not correlated)

2) If i want to include some deterministic constrains because I know that the 3 parameters i want to change are correlated, i.e. if I increase $v_1$ of $v_p$, then I want to decrease $v_2$ and $v_3$ of $v_p / 3$ (a sort of compensation)...how am I supposed to come up with a representation of the proposal ? To make it more clear: I am computing the resolution matrix for every accepted model, so I know how every parameter is correlated with the others. Suppose I want to perturb the $k^{th}$ parameter, $v_k$ of say $1km/s$. I compute the $k^{th}$ row of the resolution matrix, and I use its values to compute small perturbations that are aimed to balance the main perturbation of $v_k$.

I hope I have exposed my problem clearly, and not only in a lengthy way. Hoping that someone could help me...thanks in advance.

  • $\begingroup$ What's wrong with using vectors and density matrices/their determinants (inverse to correlation matrices) instead of numbers in the same very formula for $q(v'|v)$ you wrote? It does not look like the multivariate normal distribution is any different from the univariate one for your question unless you swept something subtle under the rug. $\endgroup$ – fedja Jul 9 '14 at 11:33
  • $\begingroup$ With constraints, all you need to do is to use conditional probabilities instead of full ones, so the "normal proposal" just gets divided by the total probability that the constraint is satisfied in the unconstrained model (I agree that the computation of this one may be ugly...). $\endgroup$ – fedja Jul 9 '14 at 11:39
  • $\begingroup$ @fedja I expanded question 2 trying to be more detailed. How can I compute the "total probability that the constraint is satisfied in the unconstrained model"? It's not so straightforward to me, sorry. $\endgroup$ – Francesco Jul 9 '14 at 12:27

Your Answer

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

Browse other questions tagged or ask your own question.