概要
環境
- python3.5
- numpy==1.13.3
- scikit-learn==0.19.1
- matplotlib==2.1.0
- scipy==1.0.0
やりたいこと
1次元、2次元のデータ分布に対して、確率密度推定を行い、その結果を以下のように図示して確認すること。
| データ分布 | 確率密度推定結果 | |
|---|---|---|
| 1D |  |
 |
| 2D |  |
 |
確率密度推定のモデル
データ標本から確率密度関数を推定する方法は、一般的に2つに大別できる。
- パラメトリックモデル: 母集団が正規分布やガンマ分布などであることを想定するモデル
- ノンパラメトリックモデル: 母集団の分布を想定しないモデル。
手法
ここでは、ノンパラメトリックの代表的な手法と言えるKDE(Kernel Density Estimation)を用いて、標本の確率密度関数を推定する。
KDEは、標本に対して複数のカーネルを適用し、それらのカーネルの重ね合わせでデータ標本の確率密度関数を推定する方法である。
KDEライブラリ
Python環境では、KDEを実装している主要なライブラリはいくつかある。
- scikit-learn
- scipy
- statsmodels
今回はscikit-learnを用いる。
1次元の場合
正規分布に従うランダムな値を生成する。
ここでは、2つの異なる正規分布を使っておく。
import numpy as np
import matplotlib.pyplot as plt
from scipy.stats.distributions import norm
np.random.seed(0)
x = np.concatenate([norm(-1, 1.).rvs(200), norm(1, 0.3).rvs(100)])
# ヒストグラムで確認する
fig, ax = plt.subplots(figsize=(17, 4))
ax.hist(x, 10, fc='gray', histtype='stepfilled', alpha=0.5, normed=True)
ax.set_xlim(x.min()-1, x.max()+1)
plt.show()
KDEを実行する。ここでは、scikit-learnにあるKDEを使う。
from sklearn.neighbors.kde import KernelDensity
X = x.reshape(-1, 1)
kde = KernelDensity(kernel='gaussian', bandwidth=0.2).fit(X)
kde.score_samples(X)
なお、kde.score_samplex()が返す値は、log-likelihood（対数尤度）であり、対数となっている。expを撮って、積分すると（ほぼ）1である。
参考までに、確認しておく。
from scipy import integrate
chk = integrate.quad(lambda x: np.exp(kde.score_samples(x)), -np.inf, np.inf)
assert np.allclose(chk, [1.0, 0.0], rtol=1e-05, atol=1e-07)
生成したランダムな値のヒストグラムとKDEの曲線を同時に図示して確認する。
x_plot = np.linspace(-5, 5, 50)[:, np.newaxis]
fig, ax = plt.subplots(figsize=(17, 4))
ax.plot(x_plot, np.exp(kde.score_samples(x_plot)), linewidth=3, alpha=0.5)
ax.hist(X, 10, fc='gray', histtype='stepfilled', alpha=0.5, normed=True)
ax.set_xlim(X.min()-1, X.max()+1)
plt.show()
2次元の場合
正規分布に従うランダムな値を2つ生成する。
ここでも、2つの異なる正規分布を使っておく。
import matplotlib.pyplot as plt
import numpy as np
from scipy.stats.distributions import norm
np.random.seed(0)
x = np.concatenate([norm(-1, 1.).rvs(200), norm(1, 0.3).rvs(100)])
y = np.concatenate([norm(-1, 1.).rvs(200), norm(1, 0.3).rvs(100)])
# 2Dヒストグラムで確認する
import matplotlib.cm as cm
fig = plt.figure(figsize=(8, 6))
ax = fig.gca()
ax.set_xlim(x.min(), x.max())
ax.set_ylim(y.min(), y.max())
H = ax.hist2d(x, y, bins=60, normed=True, cmap=cm.gray)
fig.colorbar(H[3],ax=ax)
plt.show()
KDEを実行する。
from sklearn.neighbors.kde import KernelDensity
X = np.vstack((x, y)).T
kde = KernelDensity(kernel='gaussian', bandwidth=0.2).fit(X)
np.exp(kde.score_samples(X))
ここでも、expをとって、積分すると（ほぼ）1であることを確認しておく。少しだけ時間がかかります。
from scipy import integrate
chk = integrate.dblquad(lambda _x, _y: np.exp(kde.score_samples(np.array([[_x, _y]]))), -np.inf, np.inf, gfun=lambda x: -np.inf, hfun=lambda x: np.inf)
assert np.allclose(chk, [1.0, 0.0], rtol=1e-05, atol=1e-07)
生成したランダムな値のヒストグラムとKDEの曲線を同時に図示して確認する。
f0 = np.arange(x.min(), x.max(), 0.1)
f1 = np.arange(y.min(), y.max(), 0.1)
xx, yy = np.meshgrid(f0, f1)
positions = np.vstack([xx.ravel(), yy.ravel()]).T
scores = np.exp(kde.score_samples(positions)).T
f = np.reshape(scores, xx.shape)
fig = plt.figure(figsize=(8, 6))
ax = fig.gca()
ax.set_xlim(x.min(), x.max())
ax.set_ylim(y.min(), y.max())
cset = ax.contour(xx, yy, f, colors='r')
ax.clabel(cset, inline=1, fontsize=10)
H = ax.hist2d(x, y, bins=60, normed=True, cmap=cm.gray)
fig.colorbar(H[3],ax=ax)
plt.show()
以上。