9
$\begingroup$

I have a quite flat likelihood leading Metropolis-Hastings sampler to move through the parameter space very irregularly, i.e. no convergence can be achieved no matter what the parameters of proposal distribution (in my case it is gaussian). There is no high complexity in my model - just 2 parameters, but it seems that MH cannot handle this task. So, is there any trick around this problem? Is there a sampler that would not produce Markov chains moving very far to the posterior tails?

Update of the problem:
I will try to reformulate my question giving more details. First of all I will describe the model.
I have a graphical model with two nodes. Each node is governed by an auto-Poisson model (Besag, 1974) as follows: $$p\left ( X_{j} |X_{k}=x_{k},\forall k\neq j,\Theta \right )\sim Poisson\left ( e^{\theta _{j}+\sum _{j\neq k}\theta _{kj}x_{k}} \right )$$ Or, since there just two nodes and assuming equal global intensities: $$p\left ( X_{1} |X_{2}=x_{2},\theta, \alpha \right )\sim Poisson\left ( e^{\theta+\alpha x_{2}} \right )$$ $$p\left ( X_{2} |X_{1}=x_{1},\theta, \alpha \right )\sim Poisson\left ( e^{\theta+\alpha x_{1}} \right )$$

Since it is a Markov field, the joint distribution (or likelihood of realization $ X=[x_{1},x_{2}] $) is as follows: $$ p\left ( X \right )=\frac{exp\left ( \theta \left ( x_{1}+x_{2} \right )+2 x_{1}x_{2} \alpha\right )}{Z\left ( \theta, \alpha \right )}=\frac{exp\left ( E\left ( \theta, \alpha, X \right ) \right )}{Z\left ( \theta, \alpha \right )} $$ Since I assumed flat priors for $\alpha$ and $\theta$, posterior is then proportional to $$\pi(\theta, \alpha |X)\propto \frac{exp\left ( E\left ( \theta, \alpha, X \right ) \right )}{Z\left ( \theta, \alpha \right )}$$ Since $Z(\theta, \alpha)$ in general is very hard to evaluate (lots of lots of summations) I am using auxiliary variable method due to J. Moller (2006). According to this method, first I draw a sample of data ${X}'$ by Gibbs sampler (since conditionals are just poisson distributions) then I draw a proposal from Gaussian distribution and calculate accordingly the acceptance criteria $H({X}',{\alpha}',{\theta}'|X, \alpha, \theta)$. And here I get a wild Markov chain. When I impose some boundaries within which the chain can move, the sampler seems to converge to some distribution, but once I move at least one boundary, resulting distribution also moves and always shows trancation.
I think that @Xi'an is wright - the posterior might be improper.

Improve this question
$\endgroup$
9
  • 1
    $\begingroup$ A possibility is to use a larger scale parameter in order to get larger steps. You might be interested on the R package mcmc and the command metrop as well. You will probably need an adaptive sampler. This sampler (the twalk) can be used in this kind of cases given that it is adaptive (perhaps just as a "second opinion"). It is implemented in R, C and Python. The codes can be downloaded from one of the author's webpage. $\endgroup$ – user10525 Jul 28 '12 at 19:05
  • $\begingroup$ @Procrastinator Can you elaborate more about what do you mean by "larger scale parameter"? Do mean to use larger variance parameters for proposals? $\endgroup$ – Tomas Jul 28 '12 at 19:10
  • 2
    $\begingroup$ Just let me clarify first that, if the likelihood is flat, you do not really want that your sampler do not "move very far to the posterior tail". What is desired is to sample properly from the distribution (both, tails and centre). When using a MH algorithm with Gaussian proposals, you need to choose scale parameters/covariance matrix that determine the length of the steps. These have to be chosen for 1. Sampling properly from the distribution and 2. Getting a reasonable acceptance rate. $\endgroup$ – user10525 Jul 28 '12 at 19:14
  • $\begingroup$ if you only have two parameters then numerical integration is probably a better alternative $\endgroup$ – probabilityislogic Jul 29 '12 at 1:40
  • $\begingroup$ there is something wrong with the joint likelihood expression. If you try and sum out $x_1$ you get $p(x_2|\alpha\theta)=g(x_2)\sum_{x_1=0}^{\infty}\exp(x_1[\theta+2\alpha x_2])=\infty$. so the likelihood is improper as currently written. $\endgroup$ – probabilityislogic Jul 30 '12 at 13:13
8
$\begingroup$

I find it surprising that a flat likelihood produces convergence issues: it is usually the opposite case that causes problems! The usual first check for such situations is to make sure that your posterior is proper: if not it would explain for endless excursions in the "tails". If the posterior is indeed proper, you could use fatter tail proposals like a Cauchy distribution... And an adaptive algorithm à la Roberts and Rosenthal.

If this still "does not work", I suggest considering a reparameterisation of the model, using for instance (i.e. if there is no other natural parametrisation) a logistic transform, $$ \varphi(x) = \exp(x)/\{1+\exp(x)\} $$ (with a possible scale parameter), which brings the parameter into the unit square.

Regarding the earlier answers, Gibbs sampling sounds like a more likely solution than accept-reject, which requires finding a bound and scaling the t distribution towards the posterior, which did not seem feasible for the more robust Metropolis-Hastings sampler...

Improve this answer
$\endgroup$
2
  • $\begingroup$ @Xian thanks for the feedback on the downvote. Is there actually any situation where you would favor accept-reject over MH? $\endgroup$ – gui11aume Jul 28 '12 at 23:50
  • $\begingroup$ @gui11aume: if you can produce an accept-reject algorithm with a small enough bound to ensure a reasonable acceptance rate, then accept-reject is undoubtedly preferable to Metropolis-Hastings. However, this is unlikely to happen with (a) large dimensions and/or (b) complex, possibly multimodal, targets... $\endgroup$ – Xi'an Jul 30 '12 at 10:22
2
$\begingroup$

Can you write down the distribution of your first parameter conditional on your second parameter and vice-versa? If so, Gibbs sampling would be a viable option. It's only a couple of lines of code and it can mix almost instantly in many cases.

Improve this answer
$\endgroup$
1
$\begingroup$

EDIT: See the answer of @Xi'an and the discussion after it to see the issues with the following approach.

If Metropolis-Hastings fails and your model is relatively simple, you could think of using the accept-reject algorithm with Student's $t$ distribution with a low degree of freedom (1-6) for the proposals.

If you use R, you can easily simulate a Student's $t$ with rt(). If you do not have an easy way to generate $t$ variables with your software but you can simulate a $\Gamma$, then drawing the variance of a Gaussian from a $\Gamma$ at each step and simulating a Gaussian with that variance is equivalent.

Improve this answer
$\endgroup$

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.