0
$\begingroup$

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!

$\endgroup$
2
  • $\begingroup$ 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)$? $\endgroup$ – πr8 Jan 9 at 20:35
  • $\begingroup$ @πr8 That is correct, yes ! $\endgroup$ – Kenny Abinski Jan 9 at 21:25
0
$\begingroup$

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)
curve(.5*exp(-abs(x)),add=TRUE,col="tomato")

enter image description here

The Python code posted by the OP is

  1. using a Normal mixture as proposal density and in the Metropolis-Hastings acceptance probability
  2. simulating a sum of Normal variates instead as proposed value
  3. using the sum of the logarithms of the Normal densities instead of the logarithm of the sum of the Normal densities (if the proposal is truly a mixture of two Normals, see 1. and 2. above)
$\endgroup$
5
  • $\begingroup$ Thank you for your answer! Tried to follow your code and my own, with poor results. However the code works just fine with R as you did, but Im not sure why it wont work with Python. $\endgroup$ – Kenny Abinski Jan 9 at 22:26
  • $\begingroup$ 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? $\endgroup$ – Kenny Abinski Jan 10 at 9:35
  • $\begingroup$ Or is it even possible to use sum of Normal variates as the proposal density at all? $\endgroup$ – Kenny Abinski Jan 10 at 10:56
  • $\begingroup$ 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. $\endgroup$ – Kenny Abinski Jan 11 at 11:39
  • $\begingroup$ The errors are unconnected with Markovian properties and they are listed in my answer. This is basic probability theory: the density of a sum of variates IS NOT the sum of the densities. $\endgroup$ – Xi'an Jan 11 at 11:57

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.