# Metropolis-Hastings with non-centered Proposal

I am trying to draw samples from the Laplace distribution $$\pi^* = \text{exp}(-|\theta|)$$, using Metropolis Hastings algorithm with a noncentered proposal, meaning that regular Metropolis wont work.. Now I know for a fact this is not converging, but I am really clueless as to why.

I am trying to sample it using python code, so I will use it to demonstrate my approach.

def log_target(x):
return -np.abs(x)

def eval_log_q(xp, x):
return scipy.stats.norm.logpdf(xp, x, 1) + scipy.stats.norm.logpdf(xp, x+1, 1)

def sample_q(x):
return x + npr.normal(0, 1) + npr.normal(1, 1)



where the samples are checked for acceptance as follows:

x_prop = sample_q(x)
a = np.minimum(1, np.exp(log_target(x_prop) + eval_log_q(x, x_prop) - log_target(x) - eval_log_q(x_prop, x)))

u = npr.rand()
if a > u:
#accept


I know that the chain would converge if I changed the proposal to a centered symmetric one, and thus make it ordinary metropolis sampler. But I feel that this should work as well, since this chain should be irreducible and aperiodic, thus being ergodic. Where am I wrong with this one?

Thank you!

• Just so I understand: is your proposal to set $x' = x + \xi_1 + \xi_2$, where $\xi_1 \sim N(0, 1)$, $\xi_2 \sim N(1, 1)$?
– πr8
Commented Jan 9, 2021 at 20:35
• @πr8 That is correct, yes ! Commented Jan 9, 2021 at 21:25

The proposed value need be simulated from a mixture of Normals, not as a sum of Normals. Here is an R version of the corrected code that works as a charm:

log_target<-function(x) return(-abs(x))

eval_log_q<-function(xp, x)
return(log(dnorm(xp, x, 1) + dnorm(xp, x+1, 1)))

sample_q<-function(x) return(x + (runif(1)<.5)+rnorm(1))

x=rep(rnorm(1),1e5)
for(t in 2:1e5){
x[t]=sample_q(x[t-1])
if(log(runif(1))>log_target(x[t])-log_target(x[t-1])-
eval_log_q(x[t],x[t-1])+eval_log_q(x[t-1],x[t]))
x[t]=x[t-1]}

hist(x,prob=TRUE,nclass=345)

• Ah, the sum of logarithms was indeed wrong, noticed it only now. If I wanted to keep the proposal simulation as sum of Normal variates, and thus use sum of normal ass proposal density, the code for eval_log_q would be return log(dnorm(xp, x+1, sqrt(2))) right? As the density of sum of normal variates should be $X + Y \sim N(\mu_1 + \mu_2, \sqrt{\sigma^2_1 + \sigma^2_2})$ If im correct? Commented Jan 10, 2021 at 9:35
• In R code, if my sample_q function is sample_q<-function(x) return(x + rnorm(n=1, mean=0, var=1) +rnorm(n=1, mean=1, var=1)) , would the corresponding eval_loq_q function be eval_log_q<-function(xp, x) return(log(dnorm(xp, x+1, 1)))? If not, what would the corresponding eval_log_q function be? And Im still confused as to what markov chain property does my Original Python code break to prevent converging in general. Commented Jan 11, 2021 at 11:39