I'm trying to use Metropolis-Hastings to sample from a distribution that's very close to


and I'm finding that the method is undersampling near the origin, where there's a kink in the distribution function. I've attached an example figure. Is there a smoothness requirement to Metropolis-Hastings that isn't discussed that often? If so, is there a good workaround for this or alternative algorithm?

enter image description here

Edit: Here is a minimal working example, along with the output.

def exp_func(x, ell):

    return np.exp(-np.abs(x)/ell)

ell = 2.

sigma = 1.

n_samples = 5000
sequence = np.zeros(n_samples)

theta_t = np.random.random()
idx = 0

while idx < n_samples:
    theta_star = np.random.normal(loc=theta_t, scale=sigma)
    alpha = min(exp_func(theta_star, ell)/exp_func(theta_t, ell), 1.)
    u = np.random.random()
    if u < alpha:
        theta_t = theta_star
        sequence[idx] = theta_star
        idx += 1

and what it produces:

enter image description here

  • $\begingroup$ Done. I hadn't seen anything, so I figured it might be a coding error, but as you note the fit is so good it's hard to imagine where that would be. $\endgroup$ – webb Jun 9 '20 at 21:52
  • 1
    $\begingroup$ Since the code stops there, do you repeat theta_t when u > alpha as I do in my R code (last line)? $\endgroup$ – Xi'an Jun 10 '20 at 6:03
  • $\begingroup$ Yes, the only time theta_t changes is if u < alpha. It's worth noting that I tested this with an exp(-x**4/sigma**4) type distribution and I don't get this undersampling near the origin. $\endgroup$ – webb Jun 10 '20 at 17:22
  • $\begingroup$ Yes, that's what happens. If u < alpha then I change theta_t, otherwise it stays the same, i.e. when u > alpha. $\endgroup$ – webb Jun 10 '20 at 21:49
  • 1
    $\begingroup$ I fear we have a communication problem, please check my addendum in the answer below. $\endgroup$ – Xi'an Jun 11 '20 at 6:32

Differentiability or even continuity is not a requirement for Metropolis to apply (and converge) and the fit is quite good except for the very centre, so I rather foresee a coding error.

Here is an illustration of plain Metropolis working out very well on an example from our book for the target function$$f(x) \propto \exp(-x^2/2) (\sin(6x)^2+3\cos(x)^2\sin(4x)^2+1)$$which gets easily simulated by accept-reject (Section 2.3.1).

for(t in 1:1e4){
  y=rnorm(1)#standard Normal proposal

with indeed a very good fit:

enter image description here

If one forgets to duplicate the current value under a rejection, i.e., to set $X_{t+1}=X_t$, as in

 for(t in 1:1e4){y=rnorm(1)

the fit is lost:

enter image description here

For the Laplace distribution, here is a comparison between replicating and not-replicating:

enter image description here

  • $\begingroup$ That's interesting, and yes, I saw my implementation work for other examples with smooth peaks. But your example doesn't have a discontinuity in the derivative, and mine does. Does that matter? Put in a more pathological way, what if my distribution was a step function? Is there a known problem with sampling near a discontinuity in the distribution? $\endgroup$ – webb Jun 9 '20 at 16:26

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.