適当なデータをsklearnでカーネル密度推定する

クロスバリデーション付きでどう書くのが楽かなと考えた。
私なりの答えはこんなん。ライブラリのimport設定は以下で使うものもまとめて書いちゃってる。

import numpy as np
import matplotlib.pyplot as plt
from scipy import integrate
from sklearn.model_selection import GridSearchCV
from sklearn.neighbors import KernelDensity

def estimate_kde(x, bandwidth=np.linspace(0.4, 4.0, 50), cv=5):
    # cross-validation
    grid = GridSearchCV(KernelDensity(), {'bandwidth': bandwidth}, cv=cv)
    grid.fit(x.reshape(-1, 1))
    return grid.best_estimator_

適当なサンプルデータ(ここでは混合正規分布)を作成してこれを元に実験する。
混合正規の作り方は適当なので、これよりいいやり方があったら是非知りたい。

size = 10**3
x1 = np.random.normal(0.0, 1.0, size)
x2 = np.random.normal(3.0, 0.5, size)
uniform = np.random.rand(size)
x = np.array([(x1[i] if x > 0.3 else x2[i]) for i, x in np.ndenumerate(uniform)])
# 念のための可視化
plt.hist(x)
plt.show()

カーネル密度推定を実行してみる&正しく確率密度関数になってるか、幅広い区間積分して1になるかでチェック→良さげ。
score_sample()はlogとった値返してきやがるので、expとらないとダメな点に注意、くそ仕様め。

# カーネル密度推定の実行
kde = estimate_kde(x)
# 積分してチャンと1になるか → なる
integrate.quad(lambda x: np.exp(kde.score_samples(x)), -np.inf, np.inf)

元のデータと推定した確率密度関数を重ねてみる、ここもexpとらないとダメな点に注意だ…
結果は良さげ。

# 適当なX軸の範囲で描画
x_plot = np.linspace(-5, 5, 50)[:, np.newaxis]
fig, ax = plt.subplots()
ax.plot(x_plot, np.exp(kde.score_samples(x_plot)), linewidth=3, alpha=0.5)
ax.hist(x, 30, fc='gray', histtype='stepfilled', alpha=0.3, normed=True)
ax.legend(loc='upper left')