# Simplex Random Walk

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.

## 1 Answer

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., $$x_{new} = x_{old}e^Z$$ Taking logarithms, $$\log x_{new} = \log x_{old} + Z.$$ So, actually the conditional distribution of $x_{new}$ is normal: $$\log x_{new} \mid x_{old} \sim N(\log x_{old}, h).$$ 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: $$\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}}.$$

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 $$\frac{p(x_{new})}{p(x_{old})} = \frac{e^{-x_{new}}}{e^{-x_{old}}}$$ 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.

• 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. Jun 12 '14 at 12:55