単純な 1 次元カーネル密度推定

はじめに

カーネル密度推定は、確率変数の確率密度関数を推定するために使用される統計手法です。この実験では、scikit - learn Python ライブラリを使用して、1 次元のカーネル密度推定の原理を示します。

VM のヒント

VM の起動が完了したら、左上隅をクリックしてノートブックタブに切り替え、Jupyter Notebook を使って練習しましょう。

時々、Jupyter Notebook が読み込み終了するまで数秒待つ必要があります。Jupyter Notebook の制限により、操作の検証は自動化できません。

学習中に問題に遭遇した場合は、Labby にお問い合わせください。セッション後にフィードバックを提供してください。すぐに問題を解決いたします。

ヒストグラムとカーネルの描画

まず、ヒストグラムとカーネルを描画して、両者の違いを示します。ガウスカーネル密度推定を使用して、両者の違いを示します。また、scikit - learn に用意されている他のカーネルとも比較します。

## 必要なライブラリをインポート
import numpy as np
import matplotlib.pyplot as plt
from scipy.stats import norm
from sklearn.neighbors import KernelDensity

## データを生成
np.random.seed(1)
N = 20
X = np.concatenate(
    (np.random.normal(0, 1, int(0.3 * N)), np.random.normal(5, 1, int(0.7 * N)))
)[:, np.newaxis]
X_plot = np.linspace(-5, 10, 1000)[:, np.newaxis]
bins = np.linspace(-5, 10, 10)

## グラフと軸を作成
fig, ax = plt.subplots(2, 2, sharex=True, sharey=True)
fig.subplots_adjust(hspace=0.05, wspace=0.05)

## ヒストグラム 1 を描画
ax[0, 0].hist(X[:, 0], bins=bins, fc="#AAAAFF", density=True)
ax[0, 0].text(-3.5, 0.31, "ヒストグラム")

## ヒストグラム 2 を描画
ax[0, 1].hist(X[:, 0], bins=bins + 0.75, fc="#AAAAFF", density=True)
ax[0, 1].text(-3.5, 0.31, "ヒストグラム、ビンをシフト")

## トップハット KDE を描画
kde = KernelDensity(kernel="tophat", bandwidth=0.75).fit(X)
log_dens = kde.score_samples(X_plot)
ax[1, 0].fill(X_plot[:, 0], np.exp(log_dens), fc="#AAAAFF")
ax[1, 0].text(-3.5, 0.31, "トップハットカーネル密度")

## ガウス KDE を描画
kde = KernelDensity(kernel="gaussian", bandwidth=0.75).fit(X)
log_dens = kde.score_samples(X_plot)
ax[1, 1].fill(X_plot[:, 0], np.exp(log_dens), fc="#AAAAFF")
ax[1, 1].text(-3.5, 0.31, "ガウスカーネル密度")

## データポイントを描画
for axi in ax.ravel():
    axi.plot(X[:, 0], np.full(X.shape[0], -0.01), "+k")
    axi.set_xlim(-4, 9)
    axi.set_ylim(-0.02, 0.34)

## 左列の y 軸ラベルを設定
for axi in ax[:, 0]:
    axi.set_ylabel("正規化密度")

## 下の行の x 軸ラベルを設定
for axi in ax[1, :]:
    axi.set_xlabel("x")

## グラフを表示
plt.show()

利用可能なカーネルの描画

利用可能なすべてのカーネルを描画して、それらの形状を示します。

## データを生成
X_plot = np.linspace(-6, 6, 1000)[:, None]
X_src = np.zeros((1, 1))

## グラフと軸を作成
fig, ax = plt.subplots(2, 3, sharex=True, sharey=True)
fig.subplots_adjust(left=0.05, right=0.95, hspace=0.05, wspace=0.05)

## x 軸のラベル用のフォーマット関数
def format_func(x, loc):
    if x == 0:
        return "0"
    elif x == 1:
        return "h"
    elif x == -1:
        return "-h"
    else:
        return "%ih" % x

## 利用可能なカーネルを描画
for i, kernel in enumerate(
    ["gaussian", "tophat", "epanechnikov", "exponential", "linear", "cosine"]
):
    axi = ax.ravel()[i]
    log_dens = KernelDensity(kernel=kernel).fit(X_src).score_samples(X_plot)
    axi.fill(X_plot[:, 0], np.exp(log_dens), "-k", fc="#AAAAFF")
    axi.text(-2.6, 0.95, kernel)

    axi.xaxis.set_major_formatter(plt.FuncFormatter(format_func))
    axi.xaxis.set_major_locator(plt.MultipleLocator(1))
    axi.yaxis.set_major_locator(plt.NullLocator())

    axi.set_ylim(0, 1.05)
    axi.set_xlim(-2.9, 2.9)

## 2 行目のタイトルを設定
ax[0, 1].set_title("利用可能なカーネル")

## グラフを表示
plt.show()

1 次元密度の例の描画

1 次元で 100 個のサンプルを持つ 1 次元密度の例を描画します。トップハット、ガウス、エパネチニコフの 3 種類の異なるカーネル密度推定を比較します。

## データを生成
N = 100
np.random.seed(1)
X = np.concatenate(
    (np.random.normal(0, 1, int(0.3 * N)), np.random.normal(5, 1, int(0.7 * N)))
)[:, np.newaxis]

X_plot = np.linspace(-5, 10, 1000)[:, np.newaxis]

true_dens = 0.3 * norm(0, 1).pdf(X_plot[:, 0]) + 0.7 * norm(5, 1).pdf(X_plot[:, 0])

## グラフと軸を作成
fig, ax = plt.subplots()

## 入力分布を描画
ax.fill(X_plot[:, 0], true_dens, fc="black", alpha=0.2, label="入力分布")

## 色とカーネルを設定
colors = ["navy", "cornflowerblue", "darkorange"]
kernels = ["gaussian", "tophat", "epanechnikov"]
lw = 2

## カーネル密度推定を描画
for color, kernel in zip(colors, kernels):
    kde = KernelDensity(kernel=kernel, bandwidth=0.5).fit(X)
    log_dens = kde.score_samples(X_plot)
    ax.plot(
        X_plot[:, 0],
        np.exp(log_dens),
        color=color,
        lw=lw,
        linestyle="-",
        label="kernel = '{0}'".format(kernel),
    )

ax.text(6, 0.38, "N={0} points".format(N))

## 凡例を設定してデータポイントを描画
ax.legend(loc="upper left")
ax.plot(X[:, 0], -0.005 - 0.01 * np.random.random(X.shape[0]), "+k")

## x と y の範囲を設定
ax.set_xlim(-4, 9)
ax.set_ylim(-0.02, 0.4)

## グラフを表示
plt.show()

まとめ

この実験では、scikit - learn Python ライブラリを使って、1 次元のカーネル密度推定の原理を示す方法を学びました。ヒストグラムとカーネル、利用可能なカーネル、および 1 次元密度の例を描画しました。また、異なるカーネル密度推定を比較しました。