Metropolis sampling with different proposals

I implemented a metroplis sampler for a 1D gaussian mixture, the target distribution looks like this: 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: 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
• In both cases, the convergence is slower, which means that 10⁶ iterations is presumably not sufficient to reach it. – Xi'an Jan 27 '17 at 12:43
• 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 ... – Maximal Jan 27 '17 at 14:52
• There is no mistake in your code if this is the MC core and there is no reason for "systematic bias" when running it... – Xi'an Jan 27 '17 at 15:07
• @Maximal Have you tried overlaying the true density on the plots. In my experience, sometimes aspect ratios and perspectives play with our eyes. – Greenparker Feb 2 '17 at 23:36

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. 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. 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 <- 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),
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()