Edited at

確率密度関数の推定と結果の確認


概要


環境


  • python3.5

  • numpy==1.13.3

  • scikit-learn==0.19.1

  • matplotlib==2.1.0

  • scipy==1.0.0


やりたいこと

1次元、2次元のデータ分布に対して、確率密度推定を行い、その結果を以下のように図示して確認すること。

データ分布
確率密度推定結果

1D
image.png
image.png

2D
image.png
image.png


確率密度推定のモデル

データ標本から確率密度関数を推定する方法は、一般的に2つに大別できる。


  1. パラメトリックモデル: 母集団が正規分布やガンマ分布などであることを想定するモデル

  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()

image.png

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()

image.png


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()

image.png

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()

image.png

以上。