け日記

最近はPythonでいろいろやってます

scipy.statsでカーネル密度推定 (KDE)

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度が山になっていることがわかります。

f:id:ohke:20181006212723p:plain

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

双峰性になっていることがわかります。

f:id:ohke:20181006212804p:plain

なお、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()

値が大きくなるほど、山が消え、平坦に近づくことがわかります。

f:id:ohke:20181006221618p:plain