The below python code implements the Metropolis algorithm and samples from a single variable gaussian distribution. The initial value is sampled uniformly within 5 standard deviations of the mean. Following perturbations are sampled uniformly (+/- 1 standard deviation) and added to current value. A random event is generated in range [0,1], if this value is less than the likelihood ratio of proposed/current, the movement is executed. Otherwise, the current is maintained for another iteration.

Because I'm sampling perturbations from a uniform distribution, inherently, symmetric, I'm just executing the Metropolis algorithm. I'd like to understand MH better, which makes use of (and accounts for) non-symmetric proposal distributions. A few questions:

(1) Why would we want to sample from a non-symmetric proposal distribution and can you provide a concrete example of one (which would take the place of the random.uniform(0,1) line)?

(2) Can you alter the code detailed below to change M -> MH, and make use of the proposal distribution in the answer to (1) above?

thank you!

def normal(x,mu,sigma):
    numerator = np.exp((-(x-mu)**2)/(2*sigma**2))
    denominator = sigma * np.sqrt(2*np.pi)
    return numerator/denominator
def gaussian_mcmc(hops,mu,sigma):
    states = []
    burn_in = int(hops*0.2)
    current = random.uniform(-5*sigma+mu,5*sigma+mu)
    for i in range(hops):
        movement = current + random.uniform(-1,1)
        curr_prob = normal(x=current,mu=mu,sigma=sigma)
        move_prob = normal(x=movement,mu=mu,sigma=sigma)
        acceptance = move_prob/curr_prob
        event = random.uniform(0,1)
        if acceptance > event:
            current = movement
    return states[burn_in:]

dist = gaussian_mcmc(100_000,mu=0,sigma=1)

MH Sampler

  • $\begingroup$ The Langevin algorithm is a non-symmetric example. $\endgroup$ – Xi'an Aug 15 '20 at 19:57
  • $\begingroup$ @Xi'an I just googled this! While it does answer Q1, the underlying physics/math sounds very intimidating. $\endgroup$ – mjake Aug 16 '20 at 0:14

enter image description here

If one is targeting a distribution with density $f$ over a set of $\mathbb R^k$, the Langevin algorithm (MALA) uses the gradient of the target to make the proposal: $$Y_t|X_t\sim\mathcal N_k(X_t+\eta\nabla\log f(X_t),\Omega)$$where

  • $\eta>0$ is a scale factor
  • $\nabla \log f$ denotes the gradient of $\log f$
  • $\Omega$ is a $k\times k$ covariance matrix

This proposal being assymmetric, the Metropolis-Hastings acceptance ratio is $$\dfrac{f(y_t)}{f(x_t)}\dfrac{\varphi(x_t|y_t)}{\varphi(y_t|x_t)}$$ if $\varphi(y|x)$ denotes the Normal density with mean$$x+\eta\nabla\log f(x)$$ and covariance $\Omega$.

Here is an excerpt from our book, Introducting Monte Carlo methods with R, on the matter:

One of those alternatives [to the random walk Metropolis-Hastings algorithm] is the Langevin algorithm of Roberts and Rosenthal (1998) that tries to favour moves toward higher values of the target $f$ by including a gradient in the proposal, $$ Y_t = X^{(t)} + \frac{\sigma^2}{2}\,\nabla \log f(X^{(t)}) + \sigma \epsilon_t\,, \qquad \epsilon_t\sim g(\epsilon)\,, $$ the parameter $\sigma$ being the scale factor of the proposal. When $Y_t$ is constructed this way, the Metropolis-Hastings acceptance probability is equal to $$ \rho(x,y) = \min\left\{ \dfrac{f(y)}{f(x)}\,\dfrac{g\left[(x-y)/\sigma-\sigma\,\nabla \log f(y)/2\right]} {g\left[(y-x)/\sigma-\sigma\,\nabla \log f(x)/2\right]}\,,1 \right\}\,. $$ While this scheme may remind you of stochastic gradient techniques, it differs from those for two reasons. One is that the scale $\sigma$ is fixed in the Langevin algorithm, as opposed to decreasing in the stochastic gradient method. Another is that the proposed move to $Y_t$ is not necessarily accepted for the Langevin algorithm, ensuring the stationarity of $f$ for the resulting chain.

The modification of the random walk proposal may, however, further hinder the mobility of the Markov chain by reinforcing the polarization around local modes. For instance, when the target is the posterior distribution of a Gaussian mixture model, the bimodal structure of the target is a hindrance for the implementation of the Langevin algorithm in that the local mode becomes even more attractive.

  • 1
    $\begingroup$ Thanks! My initial thoughts are "why this design?" If you take the gradient of the log of the PDF to inform how much you "nudge" X(t) when proposing Y(t), the gradient will be the smallest at the expectation and very far from it, but largest at inflection points between these extremes. I believe this will cause MALA to sample areas of both low and high prob density more often than areas of moderate prob density. (First, am I right? And second, why might this be useful?) Again, thanks!! $\endgroup$ – mjake Aug 16 '20 at 14:54

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.