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$
4
  • $\begingroup$ In both cases, the convergence is slower, which means that 10⁶ iterations is presumably not sufficient to reach it. $\endgroup$
    – Xi'an
    Commented Jan 27, 2017 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
    Commented Jan 27, 2017 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
    Commented Jan 27, 2017 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$ Commented Feb 2, 2017 at 23:36

1 Answer 1

3
$\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$
1
  • $\begingroup$ hello from 4 years later! what would be a sensible approach to make the similar simulation using Gibbs Sampler? I am trying to sample from mixture of two bivariate normals where I know the mixture distribution similar to your answer above. I am having hard time to find the appropriate conditionals. Thanks!! $\endgroup$
    – boyaronur
    Commented May 1, 2021 at 18:23

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.