Sample from a simple mixture of Gaussians

double_well.py

import numpy as np
from scipy.stats import rv_continuous
import matplotlib.pyplot as plt


class double_well_gen(rv_continuous):
    def _pdf(self, x, mu, sigma):
        sigma2 = sigma**2
        return ((np.exp(-(x - mu)**2 / (2.0 * sigma2))
                 + np.exp(-(x + mu)**2 / (2.0 * sigma2)))
                / (2.0 * np.sqrt(2.0 * np.pi * sigma2)))


double_well = double_well_gen(name='double_well').freeze(mu=2, sigma=1/2)

x = np.linspace(-5, 5, 100)
y = double_well.pdf(x)
assert np.allclose(np.trapz(y, x=x), 1.0)

samples = double_well.rvs(size=100)

plt.plot(x, y)
plt.hist(samples, density=True)
plt.xlabel('$x$')
plt.ylabel('Probability density')