nonparam_density_histogram.py

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

plt.rcParams["figure.figsize"] = (5,2)
fig, ax = plt.subplots()

pdf = lambda x: 0.7 * norm.pdf(x, 0.5, 0.1) + 0.3 * norm.pdf(x, 1, 0.1)

dx = 1e-2
xaxis = np.arange(0.25, 1.25, dx)
plt.plot(xaxis, list(map(pdf, xaxis)), label="p(x)")

N = 100
np.random.seed(0)
def sample():
    if np.random.rand() < 0.7:
        return np.random.normal(0.5, 0.1)
    return np.random.normal(1, 0.1)

samples = [sample() for _ in range(N)]
plt.plot(samples, np.zeros(len(samples)), 'x', label="samples")

delta = 0.1
bin_edges = np.arange(0.25, 1.25 + delta, delta)
ax.stem(bin_edges, [3 for _ in range(len(bin_edges))], markerfmt=' ', linefmt='k:', basefmt=' ')

histogram = np.array([sum([x >= bin_edges[i] and x < bin_edges[i+1] for x in samples]) for i in range(len(bin_edges)-1)])
normalized_histogram = (histogram / (N * delta))
plt.bar(bin_edges[:-1], normalized_histogram, align='edge', width=delta, alpha=0.3, color='g')

estimated_pdf = lambda x: histogram[np.nonzero(np.logical_and(x >= bin_edges[:-1], x < bin_edges[1:]))[0][0]]
estimated_cdf = np.cumsum([estimated_pdf(x) * delta * dx for x in xaxis])
plt.plot(xaxis, estimated_cdf, label="CDF estimate")

M = 100
def sample_from_estimate():
    toss = np.random.rand()
    return xaxis[np.nonzero(estimated_cdf >= toss)[0][0]]

samples_from_estimate = [sample_from_estimate() for _ in range(M)]
plt.plot(samples_from_estimate, np.zeros(len(samples_from_estimate)), 'x', label="samples from estimate")
plt.legend()
plt.show()