1
$\begingroup$

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.

enter image description here

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.

enter image description here

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.

$\endgroup$
2
  • 3
    $\begingroup$ 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. $\endgroup$ – jbowman Apr 12 '18 at 3:48
  • $\begingroup$ Thank you @jbowman. You've provided the resolution as well as explained why the symmetric case works and the asymmetric doesn't. $\endgroup$ – RDK Apr 12 '18 at 4:08
4
$\begingroup$

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

Histogram 1

Histogram 2

$\endgroup$
0

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.