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:
Define $g(x; \mu, \sigma) = \frac{1}{\sigma}e^{-(x-\mu)^2/2\sigma^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.
Proposal distribution is $\text{Unif}(a, b)$. Initialize $u=0$.
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.
samples_[sample_idx]
whennp.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 assignu
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$