kde.py

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

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)", color='blue')

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", color='orange')

delta = 0.1
x = 1.0
ax.stem([x - delta/2, x + delta/2], [3, 3], markerfmt=' ', linefmt='k:', basefmt=' ')

k = lambda u: 1 if abs(u) <= 0.5 else 0
sorted_samples = sorted(samples)
plt.plot(sorted_samples, [k((x_n - x) / delta) for x_n in sorted_samples], label='k((x - x_n)/delta)', color='red')

K = lambda x: sum([k((x - x_n) / delta) for x_n in samples])
p = lambda x: K(x) / (N * delta)
plt.plot(xaxis, list(map(p, xaxis)), label="p(x) estimate", color='red')

box_at_x_n = lambda x, x_n: k((x - x_n) / delta)
sum_of_N_boxes_at_x = lambda x: sum([box_at_x_n(x, x_n) for x_n in samples])
p = lambda x: sum_of_N_boxes_at_x(x) / (N * delta)
plt.plot(xaxis, list(map(p, xaxis)), label="p(x) estimate", color='red')

x = 1.0
gaussian_k = lambda u: (1/(2*math.pi)**0.5) * math.exp(-0.5 * u**2)
plt.plot(sorted_samples, [gaussian_k((x_n - x) / delta) for x_n in sorted_samples], label='k((x - x_n)/delta)', color='red')

delta = 0.01
K = lambda x: sum([gaussian_k((x - x_n) / delta) for x_n in samples])
p = lambda x: K(x) / (N * delta)
plt.plot(xaxis, list(map(p, xaxis)), label="p(x) estimate", color='red')

plt.legend()
plt.show()