scipy.statsでカーネル密度推定 (KDE) を行う方法のメモです。
カーネル密度推定は、標本データから確率密度を推定するものです。
要するにヒストグラムをなめらかにすることで、データの傾向を捉えやすくします。
2017/1/1〜2017/12/31 (365日) の東京の日別平均気温を使います。気象庁のサイトからダウンロードできます。
import numpy as np import pandas as pd import matplotlib.pyplot as plt from scipy.stats import gaussian_kde print(data) #[7.4, # 7.2, # 8.0, # 8.5, # ... # 6.2, # 3.8] plt.hist(data, bins=31, range=(0, 31)) plt.show()
5〜10度と20〜25度が山になっていることがわかります。
scipy.statsでカーネル密度推定をするためには、gaussian_kdeクラスを使います。
- インスタンス生成時にデータを渡します
- インスタンスに推定したい範囲を渡すと、密度の推定値が計算されて返されます
# gaussian_kdeのインスタンスを生成 kde = gaussian_kde(data) # 0〜31で密度推定 estimates = kde(np.linspace(0, 31, num=32)) print(estimates) # [0.00409969 0.00757112 0.0128535 0.01999006 0.02838034 0.03670379 # 0.04325381 0.04661941 0.04635048 0.04318489 0.03868037 0.03447436 # 0.03161084 0.03028665 0.03006976 0.03036525 0.03082537 0.03149862 # 0.03267868 0.03457151 0.03699876 0.03934671 0.04082662 0.04088853 # 0.03948617 0.03698289 0.03378769 0.03005188 0.02569315 0.02069462 # 0.01536906 0.01032863] # グラフで表示 plt.plot(kde(np.linspace(0, 31, num=32))) plt.xlabel('temerature [C]') plt.ylabel('kde') plt.show()
双峰性になっていることがわかります。
なお、gaussian_kdeにはbw_methodというパラメータがあります。これは各点の密度の推定に使う値の範囲 (バンド幅) を決めるためのメソッドや数値を渡すことができます。
大きな数値を渡すと、広い範囲の値を使って計算されるようになるので、全体的に鈍い形に近づいていきます。
0.05, 0.2, 0.6の3つの値で密度推定してみます。
kde_05 = gaussian_kde(data, bw_method=0.05)
kde_2 = gaussian_kde(data, bw_method=0.2)
kde_6 = gaussian_kde(data, bw_method=0.6)
plt.plot(kde_05(np.linspace(0, 31, num=32)), label='bw=0.05')
plt.plot(kde_2(np.linspace(0, 31, num=32)), label='bw=0.2')
plt.plot(kde_6(np.linspace(0, 31, num=32)), label='bw=0.6')
plt.xlabel('temerature [C]')
plt.ylabel('kde')
plt.legend()
plt.show()
値が大きくなるほど、山が消え、平坦に近づくことがわかります。
