Skip to content

Instantly share code, notes, and snippets.

@vvanirudh
Last active April 8, 2024 18:18
Show Gist options
  • Save vvanirudh/10df47c86f29698f7c9b12d8d3421ddf to your computer and use it in GitHub Desktop.
Save vvanirudh/10df47c86f29698f7c9b12d8d3421ddf to your computer and use it in GitHub Desktop.
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()
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment