parzen_window_demo2.ipynb

Kernel: Python 3 (ipykernel)

# Demonstrate a non-parametric (parzen) density estimator in 1D

# Author: Gerardo Durán Martín


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

plt.rcParams["axes.spines.right"] = False
plt.rcParams["axes.spines.top"] = False


try:
    from probml_utils import savefig
except ModuleNotFoundError:
    %pip install -qq git+https://github.com/probml/probml-utils.git
    from probml_utils import savefig


def K(u, axis=0):
    return np.all(np.abs(u) <= 1 / 2, axis=axis)


def p1(x, X, h):
    """
    KDE under a unit hypercube
    """
    N, D = X.shape
    xden, _ = x.shape

    u = ((x - X.T) / h).reshape(D, xden, N)
    ku = K(u).sum(axis=1) / (N * h**D)
    return ku


def kdeg(x, X, h, return_components=False):
    """
    KDE under a gaussian kernel
    """
    N, D = X.shape
    nden, _ = x.shape

    Xhat = X.reshape(D, 1, N)
    xhat = x.reshape(D, nden, 1)
    u = xhat - Xhat
    u = norm(u, ord=2, axis=0) ** 2 / (2 * h**2)  # (N, nden)
    px = np.exp(-u)
    if not return_components:
        px = px.sum(axis=1)

    px = px / (N * h * np.sqrt(2 * np.pi))
    return px


def main():
    data = np.array([-2.1, -1.3, -0.4, 1.9, 5.1, 6.2])[:, None]
    yvals = np.zeros_like(data)
    xv = np.linspace(-5, 10, 100)[:, None]

    fig, ax = plt.subplots(2, 2)
    # Uniform h=1
    ax[0, 0].scatter(data, yvals, marker="x", c="tab:gray")
    ax[0, 0].step(xv, p1(xv, data, 1), c="tab:blue", alpha=0.7)
    ax[0, 0].set_title("unif, h=1.0")
    # Uniform h=2
    ax[0, 1].scatter(data, yvals, marker="x", c="tab:gray")
    ax[0, 1].step(xv, p1(xv, data, 2), c="tab:blue", alpha=0.7)
    ax[0, 1].set_title("unif, h=2.0")

    # Gaussian h=1
    ax[1, 0].scatter(data, yvals, marker="x", c="tab:gray", zorder=3)
    ax[1, 0].plot(xv, kdeg(xv, data, 1), c="tab:blue", alpha=0.7, zorder=2)
    ax[1, 0].plot(xv, kdeg(xv, data, 1, True), c="tab:red", alpha=0.7, linestyle="--", zorder=1, linewidth=1)
    ax[1, 0].set_title("gauss, h=1.0")
    # Gaussian h=2
    ax[1, 1].scatter(data, yvals, marker="x", c="tab:gray", zorder=3)
    ax[1, 1].plot(xv, kdeg(xv, data, 2), c="tab:blue", alpha=0.7, zorder=2)
    ax[1, 1].plot(xv, kdeg(xv, data, 2, True), c="tab:red", alpha=0.7, linestyle="--", zorder=1, linewidth=1)
    ax[1, 1].set_title("gauss, h=2.0")

    plt.tight_layout()

    savefig("parzen_windows2.pdf", dpi=300)

    plt.show()


if __name__ == "__main__":
    main()