境界条件を伴うカーネル密度推定
カーネル密度推定と境界¶
カーネル密度推定 (kernel density estimation; KDE) はカーネル関数を用いてデータの確率密度関数を推定する手法である．Pythonの場合はscipy.stats.gaussian_kdeを使うのが主である．この記事では境界条件があるときにカーネル密度推定をどうすればいいかを扱う．問題を把握するために一様分布$U(-1,1)$からデータをサンプリングして密度推定してみよう．
import numpy as np
import matplotlib.pyplot as plt
from scipy.stats import gaussian_kde, uniform
N = 1000
np.random.seed(0)
data_uni = np.random.uniform(-1.0, 1.0, N) # 一様分布からデータ生成
kde_sp = gaussian_kde(data_uni) # KDEのinstance
x_sp= np.arange(-1.5, 1.5, 1e-2) # 密度推定する範囲
pdf_sp = kde_sp(x_sp) # 各xの点における推定した密度
rv_uniform = uniform(-1.0, 2.0)
pdf_uniform = rv_uniform.pdf(x_sp) # 真の確率密度
plt.figure(figsize=(5, 2), dpi=100)
plt.plot(x_sp, pdf_uniform, lw=2.0, label="true PDF")
plt.plot(x_sp, pdf_sp, lw=2.0, label="normal KDE")
plt.hist(data_uni, density=True, alpha=0.5, label='hist', color='0.65', edgecolor='k')
plt.xlabel("x"); plt.ylabel("Density")
plt.legend(bbox_to_anchor=(1, 1), loc='upper left')
plt.tight_layout()
結果を見るとデータ中央($x=0$)付近の密度推定は概ね上手くいっているが，境界 ($x=-1,1$) での密度は過小評価されている．このように境界条件がある場合はどうすればよいかというと，境界を越えた範囲の推定された密度を反射させる．篠本先生のWebサイト(Kernel density estimation with reflection boundary)の図が大変分かりやすい．
1次元データでの基本的な使用法¶
ライブラリを読み込む．
import kalepy as kale
kale.densityでカーネル密度推定が可能である．この際，reflect=Trueにすれば最大/最小のデータ点を境界として，反射した密度を推定する．
x_basic, pdf_basic = kale.density(data_uni, probability=True) # 通常
x_reflect, pdf_reflect = kale.density(data_uni, reflect=True, probability=True) # 反射
結果を描画する．kale.carpetという関数を使うと生データを下部に描画することができる．
plt.figure(figsize=(5, 2.5), dpi=100)
plt.plot(x_sp, pdf_uniform, lw=2.0, label="true PDF")
plt.plot(x_basic, pdf_basic, lw=2.0, label='normal KDE')
plt.plot(x_reflect, pdf_reflect, "--", lw=2.0, label='reflecting KDE')
plt.hist(data_uni, density=True, alpha=0.5, label='hist', color='0.65', edgecolor='k')
kale.carpet(data_uni, label='data')
plt.xlabel("x"); plt.ylabel("Density")
plt.legend(bbox_to_anchor=(1, 1), loc='upper left')
plt.tight_layout()
境界を明示的に設定する場合¶
明示的に境界を設定することも可能である．例えば$0<x<\infty$で定義される対数正規分布のデータを用いてみよう．
sigma = 1
data_lognorm = np.random.lognormal(sigma=sigma, size=N) # 対数正規分布からデータ生成
以下は下限を0, 上限の設定を無しにした設定である (reflect=[edge, None]の部分が該当する)．
edge = 0
x_basic2, pdf_basic2 = kale.density(data_lognorm, probability=True) # 通常
x_reflect2, pdf_reflect2 = kale.density(data_lognorm, reflect=[edge, None], probability=True) # 反射
比較のために実際の対数正規分布の確率密度(pdf)を取得しておこう．
from scipy.stats import lognorm
x_lognorm = np.arange(-1, 20, 1e-2)
rv = lognorm(sigma)
結果を描画する．
plt.figure(figsize=(4, 3), dpi=100)
plt.plot(x_lognorm, rv.pdf(x_lognorm), lw=2.0, label="true PDF")
plt.plot(x_basic2, pdf_basic2, lw=2.0, label='normal KDE')
plt.plot(x_reflect2, pdf_reflect2, "--", lw=2.0, label='reflecting KDE')
plt.hist(data_lognorm, density=True, alpha=0.5, label='hist', color='0.65', edgecolor='k')
kale.carpet(data_lognorm, label='data')
plt.xlabel("x"); plt.ylabel("Density")
plt.xlim(-1, 10)
plt.legend()
plt.tight_layout()
2次元データの場合¶
データをライブラリ内のtoyデータ生成関数を用いて生成する．
data = kale.utils._random_data_2d_03(num=1e3)
2次元の場合はkale.KDEを用いる．
kde = kale.KDE(data, reflect=[[0, None], [None, 1]])
結果を描画する．kale.dist2dを用いるとよい．
fig, axes = plt.subplots(1, 3, figsize=(8, 3), dpi=100)
axes[0].scatter(data[0], data[1], s=5, alpha=0.5)
kale.dist2d(data, ax=axes[1], cmap='Oranges') # 反射なし
kale.dist2d(kde, ax=axes[2], cmap='Greens') # 反射あり
titles = ['data', 'no reflection', 'reflection']
for ax, title in zip(axes, titles):
ax.set(xlim=[-0.5, 2.5], ylim=[-0.2, 1.2], title=title)
fig.tight_layout()
kale.cornerという各軸の密度を同時に描画する関数も用意されている．
kale.corner(kde)
plt.tight_layout()
まとめ¶
- 境界条件を伴うカーネル密度推定をする場合は境界で密度を反射させる．
- kalepyを用いると上記の実装がやりやすい．
- 前の記事 : WSI＊自己教師あり学習のメモ
- 次の記事 : 非可換ゲージ理論の頂点計算機
- 関連記事 :
