2

I am using the Metropolis-Hastings algorithm for phylogenetic inference. To do so I would like to draw the substitution matrix Q from the generalized time-reversible model.

To do so I need proposal distributions for the stationary distribution and the 6 substitution probabilities. For the stationary distributions I would use a Dirichlet, but I am not sure what to take for the substitution probabilities.

My question: Can I use a Dirichlet distribution for the substitution probabilities too?

My problem is: Dirichlet samples sum up to 1 but this is not necessarily true for the 6 parameters. The domain of the proposal distribution is only a subset of the domain of the target distribution but every single substitution probability can take all values in [0,1] (Which is everything that matters to me).

I hope you understand my question :D

Slow loris
  • 1,045
  • 6
  • 26
peri4n
  • 177
  • 6
  • You say that the substitution probabilities are not required to sum to 1, are they independent? is the sum required to be <= 1? – Nick May 30 '12 at 19:34
  • They are independent and strictly positive. Because every row in a substitution matrix should sum up to 1. There is an upper bound. So they should not be much higher then 1 . These are the only constraints. – peri4n May 31 '12 at 12:10
  • 1
    If they are independent, you could draw the substitution rates from Beta distributions. – Nick May 31 '12 at 14:44

2 Answers2

1

Having a quick look at the wikipedia page I'm not sure what the fit is here with MCMC/Metropolis-Hastings algorithm.

The models described there look like continuous time Markov chains. If that's right then you can get analytic results, and wouldn't need to resort to MCMC.

Essentially the analytic solution is to do the eigendecomposition of Q into $VLV^T$. The pmf of the process at $t$ is, then $p_i(t)=\sum_{i=1}^{4}c_i v_i e^{\lambda_i}$ where the constants, $c_i$ are obtained by solving for some initial probability vector $b$ the equation $Ec=b$.

I suspect this has arisen due to some of the terminology being the same. Unless I've misunderstood the question and the wiki reference is just a jumping-off point!

conjectures
  • 3,971
  • 19
  • 36
0

Nick is right, a Beta distribution seems like the natural choice here. In principle, you can have any proposal distribution you like. As long as the proposal distribution covers the range of the parameter distribution that you wish to sample from, then MCMC is guaranteed to converge on the right distribution. This means, however, that the Dirichlet distribution is NOT appropriate, because it does not cover all possible sets of parameters; namely it does not cover those that do not sum to one.

In practice, the proposal distribution can have a big influence on how fast you reach convergence, so you should try to pick one that proposes 'good' parameter values (near the mode of the posterior distribution) more frequently. You could adapt the Beta distribution to reflect any knowledge you may have about where the mode of the true parameter distribution may lie. If you have no prior knowledge, then you might just as well propose from a uniform distribution over the interval [0,1].

FrankD
  • 336
  • 1
  • 4