5
$\begingroup$

I implemented a metroplis sampler for a 1D gaussian mixture, the target distribution looks like this:

Target distribution

I use a 1D normal distribution as propsal, that is each candidate is sampled from a normal distribution centered around the last accepted sample with a certain standard deviation.

If I use a small variance for the proposal distribution my samples are distributed according to the target distribution although it takes a significant number of samples to achieve a good coverage. If I increase the variance of the propsal distribution my acceptance ratio rapidly drops which was expected, but also my 10^6 samples are not distributed according to the target distribution. Especially the difference between the number of samples generated for each mode of the distribution vanishes: Samples I could see why this happens in the case for a proposal distribution with small variance but I don't really understand the behaviour of the sampler in the case for $\sigma = 10$ or $100$.

Could anyone explain this behaviour to me?

for i=1:n
while 1>0
    candidate = normrnd(x_curr,sd_prop);
    p_candidate =  eval_mixture_pdf(candidate);
    if p_candidate > p_curr || rand() < (p_candidate / p_curr) 
        samples(i) = candidate;
        x_curr = candidate;
        p_curr = p_candidate;
        break;
    end
end
end
$\endgroup$
  • $\begingroup$ In both cases, the convergence is slower, which means that 10⁶ iterations is presumably not sufficient to reach it. $\endgroup$ – Xi'an Jan 27 '17 at 12:43
  • $\begingroup$ I tried 10^8 iterations with $\sigma = 10$ and it doesn't chance anything (I can not attach a third image, sorry). Of course this does not proof that it's not due to the convergence but it definitely feels more like systematic bias ... $\endgroup$ – Maximal Jan 27 '17 at 14:52
  • $\begingroup$ There is no mistake in your code if this is the MC core and there is no reason for "systematic bias" when running it... $\endgroup$ – Xi'an Jan 27 '17 at 15:07
  • $\begingroup$ @Maximal Have you tried overlaying the true density on the plots. In my experience, sometimes aspect ratios and perspectives play with our eyes. $\endgroup$ – Greenparker Feb 2 '17 at 23:36
2
$\begingroup$

I tried to replicate your results and here is what happened. I think I guessed the mixture distribution almost correctly: $$.3* N(0, 2) + .7*N(10, 2)\,. $$

I implemented the scenarios that you demonstrate. Here are my various settings:

  • Monte Carlo samples of $N = 10^4$ and $N = 10^6$
  • Four $\sigma$s: .1, 1, 10, 100
  • Two starting values: 0 and 10

Here are the results from $N = 10^4$. Black line is the true density. Things are a little crazy here, but I don't see what you see.

enter image description here

Now here is the more stable $N = 10^6$. Again, I don't see what you see, and rather at $\sigma = .1$ things are a little weird, rest are all ok.

enter image description here

So I am not sure why you are seeing a weird behavior. It may be a starting value issue, but I doubt it. Your algorithm looks correct, so there might be a problem with how you're storing the Markov chain (maybe?). Here is the R code for these plots (that is the only language I know).

set.seed(100)
## Evaluates the mixture density
eval_dense <- function(x)
{
  f_x <- .3*dnorm(x, mean = 0, sd = 2) + .7*dnorm(x, mean = 10, sd = 2)
  return(f_x)
}

# Samples and plots 
plot_MC <- function(sd, N = 1e6, start = 0)
{
  out <- numeric(length = N)
  accept <- 0
  out[1] <- start

  for(i in 2:N)
  {
    cand <- rnorm(1, mean = out[i-1], sd = sd)
    # Calculate ratio
    ratio <- eval_dense(cand)/eval_dense(out[i-1])
    if(runif(1) < ratio)
    {
      out[i] <- cand
      accept <- accept+1
    }else{
      out[i] <- out[i-1]
    }
  }

  x <- seq(-10, 20, length = 10000)
  truth <- eval_dense(x)
  plot(x, truth, type = 'l',
       main = paste("N =", N, ", Sigma = ", sd,", Accept = ", accept/N,  ", Start = ", out[1]),
       ylab = "density (red is MCMC)")
  lines(density(out), col = "red")
}

png("mix_mh_4.png", height = 1000, width = 800)
par(mfrow = c(4,2))
plot_MC(sd = .1, start = 0, N = 1e4)
plot_MC(sd = .1, start = 10, N = 1e4)
plot_MC(sd = 1, start = 0, N = 1e4)
plot_MC(sd = 1, start = 10, N = 1e4)
plot_MC(sd = 10, start = 0, N = 1e4)
plot_MC(sd = 10, start = 10, N = 1e4)
plot_MC(sd = 100, start = 0, N = 1e4)
plot_MC(sd = 100, start = 10, N = 1e4)
dev.off()

png("mix_mh_6.png", height = 1000, width = 800)
par(mfrow = c(4,2))
plot_MC(sd = .1, start = 0)
plot_MC(sd = .1, start = 10)
plot_MC(sd = 1, start = 0)
plot_MC(sd = 1, start = 10)
plot_MC(sd = 10, start = 0)
plot_MC(sd = 10, start = 10)
plot_MC(sd = 100, start = 0)
plot_MC(sd = 100, start = 10)
dev.off()
$\endgroup$

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.