Metropolis -> Metropolis-Hastings for asymmetric proposal distributions?

Question

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):
        states.append(current)
        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)
plt.hist(dist,normed=1,bins=20) 
plt.plot(lines,normal_curve)

Answer 1

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.