# Convergence issue in simple 1D Metropolis algorithm

I want to write a Metropolis sampler to sample independent rvs $x$ from the mixture model $X \sim \frac{1}{2}\big[\mathscr{N}(\mu_1, \sigma_1) + \mathscr{N}(\mu_2, \sigma_2)\big]$.

My algorithm is as follows:

1. Define $g(x; \mu, \sigma) = \frac{1}{\sigma}e^{-(x-\mu)^2/2\sigma^2}$

2. Define $h(x) = g(x; \mu_1, \sigma_1) + g(x; \mu_2, \sigma_2)$. Note that this is proportional to our desired sample distribution.

3. Proposal distribution is $\text{Unif}(a, b)$. Initialize $u=0$.

4. For each sample

a. Draw candidate from proposal distribution

b. Compute $p = \min\left(1, \frac{h(v)}{h(u)}\right)$

c. If $r\sim\text{Unif(0,1)} < p$ add $v$ to samples, else reject.

d. Set $u \leftarrow v$

Below is the code:

import numpy as np
from scipy.stats import norm
import matplotlib.pyplot as pl

def h(x, mu1, sig1, mu2, sig2):
v1 = norm.pdf(x, loc=mu1, scale=sig1)
v2 = norm.pdf(x, loc=mu2, scale=sig2)
return 0.5 * (v1 + v2)

# Convenience function, does not participate in MCMC simulation.
def compute_analytical(mu1, sig1, mu2, sig2):
x = np.linspace(-2, 8, 1000)
y = h(x, mu1, sig1, mu2, sig2)
return x, y

def compute_mcmc(mu1, sig1, mu2, sig2):
args = mu1, sig1, mu2, sig2
num_samples = 10000
low, high = -10, 10
# Draw samples
samples_ = np.zeros(num_samples)
sample_idx = 0
u = 0.0
while sample_idx < num_samples:
v = np.random.uniform(low, high)
pval = min(1.0, h(v, *args) / h(u, *args))
if np.random.random() < pval:
samples_[sample_idx] = v
u = v
sample_idx += 1

return samples_

if __name__ == '__main__':
mu1, sig1 = 1.0, 1.0
mu2, sig2 = 4.0, 0.25

samples = compute_mcmc(mu1, sig1, mu2, sig2)
x_theory, y_theory = compute_analytical(mu1, sig1, mu2, sig2)

fig, axes = pl.subplots(1, 1)
axes.hist(samples, bins=25, histtype='stepfilled', normed=True, color='orange', alpha=0.5)
axes.plot(x_theory, y_theory, 'b', lw=1)
pl.show()

The algorithm seems to work when $\sigma_1 = \sigma_2$ (as judged by equal peak heights in the sample histogram). The graph below shows the histogram of the samples and the theoretical distribution.

But the algorithm seems to fail when $\sigma_1 \neq \sigma_2$ (relative peak heights are different from theoretical distribution). I've tried to generate $10\times$ more samples, but did not observe any movement toward correct peak ratios.

It is easy to see that our desired distribution is the stationary distribution of the modified Markov chain. So I'm not seeing why my implementation is failing for mixture of unequal variances.

• What gets assigned to samples_[sample_idx] when np.random.random() >= pval matters too; in your case, it doesn't look like anything gets assigned, leaving a default of 0, which doesn't seem to match the histograms. You should assign u in that case. The problems this causes will be symmetric in the first example but will reduce the height of the rightmost peak in the second; once you are near 4, most of your generated acceptance probabilities will be low, but you aren't reflecting that by capturing repeated values near 4 from failures to accept. Apr 12, 2018 at 3:48
• Thank you @jbowman. You've provided the resolution as well as explained why the symmetric case works and the asymmetric doesn't.
– RDK
Apr 12, 2018 at 4:08

Even the first plot is slightly off: the 'peaks' are being undersampled. There's a piece missing in the algorithm you used. It should be:

1. Define $h(x| \mu_1, \sigma_1, \mu_2, \sigma_2) = \frac{1}{2}\Big(\frac{1}{\sigma_1}\text{exp}\big(-\frac{(x-\mu_1)^2}{2 \sigma_1^2}\big) + \frac{1}{\sigma_2}\text{exp}\big(-\frac{(x-\mu_2)^2}{2 \sigma_2^2}\big)\Big)$
2. Use proposal distribution $u \sim \text{Unif}(-10, 10)$, initialize $u_0 = 0$
3. While $\text{i} < \text{sample_size}$:
• Draw candidate, $\text{v}$, from the proposal distribution
• Compute $\text{p} = \text{min}\big(1, \frac{h(v)}{h(u)}\big)$
• If $\text{r} \sim \text{Unif}(0,1) < p$, then $\text{sample}_i = \text{v}; \text{set u} \rightarrow \text{v}$
• Else $\text{sample}_i = \text{u}$
• Increment $\text{i}$

Note the addition of the "Else $\text{sample}_i = \text{u}$" piece. If you only add samples where $\text{r} < \text{p}$ is satisfied, then you will not draw enough samples from the highest areas of $h(x)$. If you change your compute_mcmc() function to this:

def compute_mcmc(mu1, sig1, mu2, sig2):
args = mu1, sig1, mu2, sig2
num_samples = 10000
low, high = -10, 10
# Draw samples
samples_ = np.zeros(num_samples)
sample_idx = 0
u = 0.0
while sample_idx < num_samples:
v = np.random.uniform(low, high)
pval = min(1.0, h(v, *args) / h(u, *args))
if np.random.random() < pval:
samples_[sample_idx] = v
u = v
sample_idx += 1
else:
samples_[sample_idx] = u
sample_idx += 1

return samples_

then the plots will verify the correct distribution (I increased the sample size to 30k and the number of bins in the histogram to 40 to show more detail).