# Metropolis-Hastings undersampling near kink in distribution function

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

$$\exp(-|x|/\ell)$$

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?

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:

• 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. – webb Jun 9 '20 at 21:52
• Since the code stops there, do you repeat theta_t when u > alpha as I do in my R code (last line)? – Xi'an Jun 10 '20 at 6:03
• 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. – webb Jun 10 '20 at 17:22
• Yes, that's what happens. If u < alpha then I change theta_t, otherwise it stays the same, i.e. when u > alpha. – webb Jun 10 '20 at 21:49
• I fear we have a communication problem, please check my addendum in the answer below. – 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).

x=c()
for(t in 1:1e4){
y=rnorm(1)#standard Normal proposal
x=c(x,ifelse(runif(1)<f(y)*dnorm(z)/f(z)/dnorm(y),z<-y,z))}


with indeed a very good fit:

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)
if(runif(1)<f(y)*dnorm(z)/f(z)/dnorm(y))x=c(x,z<-y)}


the fit is lost:

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

• 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? – webb Jun 9 '20 at 16:26