1
$\begingroup$

This link describes how to perform a random walk on the simplex using the Metropolis-Hastings algorithm:

http://en.wikipedia.org/wiki/User:Skinnerd/Simplex_Point_Picking

The description says:

"The stationary distribution of its samples is the unit-exponential distribution."

As well as

"This procedure effectively samples x_new from a gamma random variable with mean of x_old and standard deviation of h*x_old."

Finally the actual update step is:

"x_new <- x_old * exp( Random_Normal(0,h) )"

I'm trying to understand the algorithm, specifically how it calculates the hastings ratio given the general description of the algorithm here http://en.wikipedia.org/wiki/Metropolis%E2%80%93Hastings_algorithm.

My question is what is the proposal distribution g(x_old -> x_new)?

Is the proposal distribution actually g(x_old -> x_new) = x_old * exp( Random_Normal(0,h) ), which would make the hastings ratio g(x_new -> x_old) / g(x_old -> x_new) = x_new / x_old?

Thanks for any help.

$\endgroup$
6
$\begingroup$

The linked algorithm seems to actually be slightly erroneous (I explain the error in the end of the answer), however, the Hastings ratio you ask about is computed correctly

Goal of the algorithm

The algorithm described is a Metropolis-Hastings sampling on one-dimensional $x$ targeting the unit-exponential distribution. To use this to obtain the samples on the simplex, one needs to normalize a vector of samples from $n$ different unit-exponential $x$:s. The function next_point performs a single Metropolis-Hastings step in these one-dimensional distributions

Proposal distribution and Hastings ratio

The proposal is made in the following line:

x_new <- x_old * exp( Random_Normal(0,h) )

So, $x_{new}$ is drawn by multiplying $x_{old}$ by a $e^Z$ where $Z\sim N(0,h)$. i.e., \begin{equation} x_{new} = x_{old}e^Z \end{equation} Taking logarithms, \begin{equation} \log x_{new} = \log x_{old} + Z. \end{equation} So, actually the conditional distribution of $x_{new}$ is normal: \begin{equation} \log x_{new} \mid x_{old} \sim N(\log x_{old}, h). \end{equation} From the last expression we see that $q(x_{new}\mid x_{old})$ is a log-normal distribution with parameters $(\log x_{new},h)$. Then, to evaluate the Hastings ratio, we just plug in the log-normal densities: \begin{equation} \left. \frac{q(x_{old} \mid x_{new})}{q(x_{new} \mid x_{old})} = \frac{1}{x_{old}\sqrt{2\pi}h}e^{-\frac{-(\log x_{old} - \log x_{new})^2}{2h^2}} \middle/ \frac{1}{x_{new}\sqrt{2\pi}h}e^{-\frac{-(\log x_{new} - \log x_{old})^2}{2h^2}}\right. = \frac{x_{new}}{x_{old}}. \end{equation}

This part corresponds to the code line:

hastings_ratio <- ( x_new / x_old )

Metropolis ratio (ratio of target densities)

Finally, we need to compute the ratio between the densities of the target distribution \begin{equation} \frac{p(x_{new})}{p(x_{old})} = \frac{e^{-x_{new}}}{e^{-x_{old}}} \end{equation} which corresponds to the code line:

metropolis_ratio <- exp(-x_new) / exp(-x_old)

The error in the algorithm

The acceptance ratio is calculated correctly. However, in Metropolis-Hastings, if a proposal is rejected, the current value should be repeated in the chain. The next_point function as written in the link will retry proposing new values until a proposal is accepted. Instead, it should simply return $x_{old}$ if the proposed value is rejected.

$\endgroup$
  • $\begingroup$ Thanks, your answer has completely cleared up my confusion. I have accepted it but cannot upvote because I have less than 15 reputation on this forum. $\endgroup$ – Richy Jun 12 '14 at 12:55

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.