Convergence issue in simple 1D Metropolis algorithm

Question

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.

Answer 1

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