カーネル密度推定とは
カーネル密度推定(Kernel Density Estimation; KDE)とは、サンプルデータから確率密度関数(Probability Density Function; PDF)を推定するノンパラメトリックな方法である。ノンパラメトリックとは、明確な関数形を仮定せずにデータ分布を柔軟に推定する方法である。
直感的なイメージ
サンプルデータが以下の様にあったとする:
x = [2.0, 3.5, 4.0, 5.5, 6.0]
KDEでは、この各データ点に対して、小さな山(ガウス分布など)を配置する:
x = 2.0 : Gaussian(2.0, \sigma_{2.0}^{2.0})
x = 3.5 : Gaussian(3.5, \sigma_{3.5}^{2.0})
・・・
これらを全部重ねると、滑らかな確率密度関数が得られる。この小さな山のことをカーネル関数と呼んでいる
KDEの数学的表現
1次元の場合、カーネル密度関数は以下の様に表現される:
\hat{p(x)} = \frac{1.0}{nh}\sum_{i=1}^{n}K(\frac{x-x_i}{h})
h:バンド幅
x_i:サンプルデータ
K(u):カーネル関数
n:サンプル数
KDEのカーネル関数とは
K(u)=\frac{1}{\sqrt{2\pi}}e^{-\frac{u^{2}}{2}}
先ほどは、小さな山のことだと説明したが、この山をガンマ分布で表現するならGaussianカーネルと呼ばれる(上式)。これ以外にもRectangularカーネル、Triangleカーネル、Epanechnikovカーネルがある。典型的なカーネル関数はGaussianカーネルだと思われる。
バンド幅
バンド幅hは、推定される密度関数の形状を決定するパラメータである。バンド幅が小さい(hが小さい)と密度関数は鋭くなり、局所的な変動に敏感になる。一方で、バンド幅が大きい(hが大きい)と、密度関数は滑らかになり、大きな構造のみが捉えられる。実際に描画してみると下図のようになる。
上図から、バンド幅の選択は、KDEから得られる推定値に強く影響するため、非常に重要になる。バンド幅の決め方としては、ヒューリスティックな手法や交差検証、プラグイン法など様々ある。詳細は「D.M. Bashtannyk and R.J. Hyndman, “Bandwidth selection for kernel conditional density estimation”, Computational Statistics & Data Analysis, Vol. 36, pp. 279-298, 2001.」を参照してもらいたい。
ヒューリスティックな手法の1つであるSilvermanの方法について簡単に述べる。サンプル数と標準偏差だけで、以下の様にして求められる。しかし、これは経験則であり、精度を求めているのなら交差検証をされるべきである。
h=\frac{1.06\hat{\sigma}}{n^{\frac{1}{5}}}
上のSilvermanの式は1次元の場合である。多次元の場合は以下の様に表される。
h = \frac{4}{(d+2)n}^{\frac{1}{d+4}}\hat{\sigma}
ここで、dは次元数を意味する。また、ここでの標準偏差は分散共分散行列の対角要素の平方根である。即ち、次元数だけhが用意される。
2次元のカーネル密度推定
上記のSilvermanの式を用いて、2次元のカーネル密度推定の結果を下図に示している。x方向のバンド幅は0.494なのに対して、y方向のバンド幅は0.097である。つまり、x方向にはやや広く、y方向には非常に狭くスムージングしていることを意味する。
参考までに、上記の図を描画するのに作成したコードを以下に示している。
# --- request module
import numpy as np
import matplotlib.pyplot as plt
from scipy.stats import gaussian_kde
# --- Generate synthetic 2D data
np.random.seed(0)
x_data = np.random.normal(loc=35, scale=5, size=1000)
y_data = np.random.normal(loc=3, scale=1, size=1000)
data_2d = np.vstack([x_data, y_data])
# --- KDE with Silverman's rule of thumb
kde = gaussian_kde(data_2d, bw_method='silverman')
bandwidth_factor = kde.factor
covariance_matrix = kde.covariance
bandwidth_matrix = kde.factor ** 2 * kde.covariance
bandwidth_vector = np.sqrt(np.diag(bandwidth_matrix))
print("\nBandwidth factor:\n" , bandwidth_factor)
print("\nCovariance matrix:\n", covariance_matrix)
print("\nBandwidth matrix:\n" , bandwidth_matrix)
print("\nBandwidth vector:\n" , bandwidth_vector)
# --- Create grid to evaluate KDE
xgrid = np.linspace(15, 55, 100)
ygrid = np.linspace(-1, 7, 100)
X, Y = np.meshgrid(xgrid, ygrid)
positions = np.vstack([X.ravel(), Y.ravel()])
Z = kde(positions).reshape(X.shape)
# --- Plot KDE heatmap
plt.figure(figsize=(9, 8))
plt.contourf(X, Y, Z, levels=50, cmap='viridis')
plt.colorbar(label='Density')
plt.scatter(x_data, y_data, s=5, c='white', alpha=0.3, label='Samples')
plt.title(f"2D KDE using Silverman's Rule,\n Bandwidth={bandwidth_vector}")
plt.xlabel("X")
plt.ylabel("Y")
plt.legend()
plt.grid(True)
plt.tight_layout()
plt.show()
Scipyは非常に便利である。モジュールの詳細は"https://docs.scipy.org/doc/scipy/reference/generated/scipy.stats.gaussian_kde.html#id2 "を参考にしてもらいたい。
まとめ
ノンパラメトリックな確率密度の推定手法としてカーネル密度推定を説明した。詳細な原理が分かっていなくても、使えてしまうのが便利である。しかし、それは危険であると同時に統計の面白みが失われる。以下の文献は、密度推定を学ぶのに良書だと思われる。全てを読んだわけではないが、、「B.W. Silverman, “Density Estimation for Statistics and Data Analysis”, Vol. 26, Monographs on Statistics and Applied Probability, Chapman and Hall, London, 1986. ( https://doi.org/10.1201/9781315140919 )」