分析ノート

今回は statsmodelsを使ったカーネル密度推定を紹介します。

ドキュメントはこちらです: statsmodels.nonparametric.kde.KDEUnivariate

statsmodelsのいつもの作法で、インスタンスを作る時にデータを渡し、fitメソッドを呼び出せば完了です。
(SciPyではfitは必要なかったので注意)

fit する時に、カーネル関数の種類やバンド幅の決定方法を指定できます。
詳しくは fit メソッドのドキュメントを確認してください。

カーネル関数が 7種類指定できたり、バンド幅の指定方法が少し違ったりします。

ではやっていきましょう。 比較のため、 SciPyで推定した結果も載せていきます。
(SciPyのバンド幅指定のルールがおかしいと言う話もあるので。参考記事はこちら。)
まずはデータ生成です。

import numpy as np
import statsmodels.api  as sm
import matplotlib.pyplot as plt
from scipy.stats import norm, gaussian_kde

# 推定したい分布の真の確率密度関数
def p_pdf(x):
    return (norm.pdf(x, loc=-4, scale=2)+norm.pdf(x, loc=2, scale=1))/2


# その分布からのサンプリング
def p_rvs(size=1):
    result = []
    for _ in range(size):
        if norm.rvs() < 0:
            result.append(norm.rvs(loc=-4, scale=2))
        else:
            result.append(norm.rvs(loc=2, scale=1))
    return np.array(result)


# グラフ描写のため真の確率密度関数の値を取得しておく
xticks = np.linspace(-10, 5, 151)
pdf_values = p_pdf(xticks)

# 100件のデータをサンプリング
data = p_rvs(100)

さて、ここから SciPyと statsmodelsでそれぞれのScottのルールとSilvermanのルールでカーネル密度推定して結果を可視化します。

# SciPyのカーネル密度推定
kde_scipy_scott = gaussian_kde(data, bw_method="scott")
kde_scipy_silverman = gaussian_kde(data, bw_method="silverman")

# statusmodels のカーネル密度推定
kde_sm_scott = sm.nonparametric.KDEUnivariate(data)
kde_sm_scott.fit(kernel="gau", bw="scott")
kde_sm_silverman = sm.nonparametric.KDEUnivariate(data)
kde_sm_silverman.fit(kernel="gau", bw="silverman")

fig = plt.figure(figsize=(7, 7), facecolor="w")
ax = fig.add_subplot(1, 1, 1, title="カーネル密度推定の結果")
ax.plot(xticks, pdf_values, label="真の分布")
ax.plot(xticks, kde_scipy_scott.evaluate(xticks), label="SciPy - Scott")
ax.plot(xticks, kde_scipy_silverman.evaluate(xticks), alpha=0.5, label="SciPy - Silverman")
ax.plot(xticks, kde_sm_scott.evaluate(xticks), linestyle=":", label="statsmodels - Scott")
ax.plot(xticks, kde_sm_silverman.evaluate(xticks), label="statsmodels - Silverman")
ax.hist(data, bins=20, alpha=0.1, density=True, label="標本")
ax.legend()
plt.show()

出力がこちらです。

statsmodels で Silverman のルールを使ったのが一番良さそうですね。
ただ、これは常にそうなるわkではなく、元の分布の形や標本サイズによって何が最適かは変わってきます。

SciPyのSilvermanのルールと、statsmodelsのScottのルールが重なってしまうのは、先日の記事で書いた通り、SciPyの不備と思われます。

さて、こうしてみると、 statsmodels が一番良さそうに見えますが、当然これにも欠点があります。
SciPyやscikit-learnのカーネル密度推定では多次元のデータも扱えるのですが、
statsmodelsでは1次元のデータしか扱えません。

カーネル密度推定に関しては SciPy, scikit-learn statusmodels　それぞれにメリットデメリットあるので、
目的に応じて使い分けていくのが良さそうです。
(これにそこまで高い精度が求められることも少ないので一番手軽に書ける SciPyで十分なことが多いと思います。)

今回もカーネル密度推定です。前回の記事の予告通り、今回はscikit-learnを使います。

ドキュメントはこちらです: sklearn.neighbors.KernelDensity

scikit-learnの他のモデルと同様に、import してインスタンス作って、fitしたら学習完了、となるのですが、その後の推定に注意することがあります。
version 0.23.1 時点では、他のモデルにある、predict やtransform がありません。
これは代わりに、score_samples を使います。そして戻り値は、
Evaluate the log density model on the data.
とある様に、その点での密度の”対数”です。なので、SciPyの時みたいに密度関数を得たかったら指数関数で変換する必要があります。

さて、scikit-learn でカーネル密度推定をやるメリットですが、ガウスカーネルだけでなく、全部で6種類のカーネル関数が使えます。

1点だけのデータにたいしてカーネル密度推定を実施すると、そのカーネルの形がそのまま残るのでそれを利用してカーネル関数を可視化しました。

import numpy as np
import matplotlib.pyplot as plt
from sklearn.neighbors import KernelDensity

kernel_names = [
        "gaussian",
        "tophat",
        "epanechnikov",
        "exponential",
        "linear",
        "cosine",
    ]

# x軸のメモリようの配列
xticks = np.linspace(-3, 3, 61)

fig = plt.figure(figsize=(10, 7), facecolor="w")
for i, kernel in enumerate(kernel_names, start=1):
    # 指定したカーネルでモデル作成
    kde_model = KernelDensity(bandwidth=1, kernel=kernel)
    # 1点だけのデータで学習
    kde_model.fit([[0]])
    ax = fig.add_subplot(2, 3, i, title=kernel)
    ax.set_ylim([0, 1.1])
    ax.plot(xticks, np.exp(kde_model.score_samples(xticks.reshape(-1, 1))))

plt.show()

結果がこちら。 それぞれ個性ありますね。

さて、カーネル密度推定はscikit-learnで行う場合もバンド幅の決定の問題があります。
そして、SciPyの様に自動的には推定してくれません。

そこで、scikit-learnでは、他のモデルと同じ様に、最適なパラメーターを探索することになります。
Score というメソッドで、与えられたデータに対する対数尤度の合計値が算出できるので、これを使ってグリッドサーチすると良さそうです。
他のどんなモデルもそうなのですが、訓練データと評価データはかならず分ける必要があります。
特にカーネル密度推定では、バンド幅が狭ければ狭いほど訓練データに対するスコアは高くなり、極端に過学習します。
今回は久々にグリッドサーチでもやりましょうか。
また、正解の分布はベータ分布を使います。 (正規分布だとガウスカーネルが有利すぎる気がしたので)

先ほどの ６種類のカーネル + SciPyでやってみました。

from scipy.stats import beta
from sklearn.model_selection import GridSearchCV

# グリッドサーチでバンド幅決めてベストのモデルを返す関数
def search_bandwidth(data, kernel="gaussian", cv=5):
    bandwidth = np.linspace(0.1, 2, 190)
    gs_model = GridSearchCV(
                KernelDensity(kernel=kernel),
                {'bandwidth': bandwidth},
                cv=cv
            )
    gs_model.fit(data.reshape(-1, 1))
    return gs_model.best_estimator_

# 正解の分布
frozen = beta(a=3, b=2, loc=-1, scale=2)
# データ生成
data = frozen.rvs(200)

# x軸のメモリようの配列
xticks = np.linspace(-1.1, 1.1, 221)


fig = plt.figure(figsize=(10, 11), facecolor="w")
for i, kernel in enumerate(kernel_names, start=1):
    # 指定したカーネルでモデル作成
    kde_model = search_bandwidth(data, kernel=kernel)
    bw_str = str(kde_model.bandwidth.round(2))
    ax = fig.add_subplot(3, 3, i, title=kernel + " h=" + bw_str)
    # ax.set_ylim([0, 1.1])
    ax.plot(xticks, np.exp(kde_model.score_samples(xticks.reshape(-1, 1))))
    ax.plot(xticks, frozen.pdf(xticks), linestyle="--")

# おまけ。scipy
kde = gaussian_kde(data)
ax = fig.add_subplot(3, 3, 7, title="SciPi")
ax.plot(xticks,  kde.evaluate(xticks))
ax.plot(xticks, frozen.pdf(xticks), linestyle="--")

plt.show()

結果がこちら。 一番手軽に使える SciPyが頑張ってますね。
全データを学習に使えるのも大きいかもしれません。

僕の環境の SciPyは version 1.4.1 なのでこの記事ではそれを想定して読んでください。(将来のバージョンでは修正されるかも。)

さて、カーネル密度推定においてバンド幅(パラメーターh)を決めるのに使われる、
Scottのルールと、Silvermanのルールについて紹介しようと思っていたのですが、
SciPyの実装を精査していく中でちょっと挙動がおかしいことに気付きました。

それぞれのルールの原典がどっちも無料では見れないのでそれを確認できないのですが、
SciPyで実装されているSilvermanのルールが、どうやら他で言うScottのルールのようです。
SciPyで実装されているScottのルールは、これ一体何なんでしょう。
他で言われているSilvermanのルールはもっと狭いバンド幅で推定するはずです。

Scottのルールの原典:
D.W. Scott, “Multivariate Density Estimation: Theory, Practice, and Visualization”, John Wiley & Sons, New York, Chicester, 1992.
Silvermanのルールの原典:
B.W. Silverman, “Density Estimation for Statistics and Data Analysis”, Vol. 26, Monographs on Statistics and Applied Probability, Chapman and Hall, London, 1986.

本当は正しい(と言うか、この二つのルールも経験則的なルールであって、それを正しいと呼ぶかは疑問符が付きますが)それぞれの発案者が提案したルールを紹介したいのですが、
普段の利用を考えると、scipy.stats.gaussian_kde で、bw_method　に指定した文字列によって挙動がどう変わるのかを正確に理解している方が役に立つと思います。
そのため、この記事のタイトル通り、発案者が何を提唱したかは一旦脇に置いといて、文字列として、”scott”と”silverman”を指定したらプログラムがどう計算するかを紹介します。

ドキュメントはこちらです: scipy.stats.gaussian_kde
GitHubのコードもかなり参考になりました。

さて、推定したバンド幅はその２乗の値が covariance という属性に格納されています。(１次元データの場合。)
これは、 学習したデータの不変分散に、 factor という属性の”2乗を”掛けることで計算される実装になっています。
バンド幅hで考えると、不変分散の平方根に factor を掛けたものになりますね。

で、その factor の計算方法が、指定した “scott”と”silverman”で変わります。
データの件数を$n$、データの次元を$d$とすると、
“scott”を指定したときは、
$$
factor = n^{\frac{-1}{d+4}}
$$
です。1次元の場合は$d=1$なので$n^{-0.2}$ですね。

“silverman”を指定したときは、
$$
factor = \left(\frac{d+2}{4}*n\right)^{\frac{-1}{d+4}}
$$
になります。
$d=1$の時は、$(3/4)^{-0.2}=1.05922\cdots$なので、だいたい、$1.06n^{-0.2}$くらいになります。
他のサイトの情報を見ていると、 0.9とか1.06とかの定数が出てくるのですが、そのうち1.06がどこから登場するのかはこれで分かりますね。
ちなみに statsmodelsのカーネル密度推定では 1.059がハードコーディングされています。

で、この 1.06は本当は Scott のルールで使われる定数のはずなのですが、なぜかSciPyでは”silverman”を指定したときに登場します。
1.06の代わりに0.9くらいの値を使うのがSilvermanのルールのはずなのですが。

ただ、以上で、SciPyがどの様にバンド幅を決定しているかはわかったので、その通りの挙動をしていることをコードで確認しておきましょう。
とりあえずデータを準備します。

from scipy.stats import gaussian_kde
from scipy.stats import norm
import numpy as np

# 50件のデータをサンプリング
data = norm(loc=2, scale=2).rvs(50)

print("標本の不偏分散: ", np.var(data, ddof=1))
# 標本の不偏分散:  4.2218467749903645

まずは “scott” の確認です。(初期値なので、コード中に”scott”の文字列は登場しません。)

# Scottのルールでカーネル密度推定
kde = gaussian_kde(data)

# バンド幅 hの2乗を確認
print("バンド幅 h^2: ", kde.covariance[0, 0])
# バンド幅 h^2:  0.8829059945819668

# Scottの因子が、データ件数の -0.2乗に等しいことを確認
print(kde.scotts_factor())
# 0.45730505192732634
print(len(data) ** (-0.2))
# 0.45730505192732634

# データの普遍分散に Scottの因子の2乗を掛けてみると、先ほどのバンド幅の２乗に等しい
print((kde.scotts_factor()**2) * np.var(data, ddof=1))
# 0.8829059945819668

以上の通り、ドキュメント通りの実装になっていました。

次が bw_method=”silverman” の場合です。ごくわずか謎の端数が出てますが、こちらもドキュメント通りの実装であることが確認できます。

# カーネル密度推定 (silvermanのルール)
kde = gaussian_kde(data, bw_method="silverman")

# バンド幅 hの2乗を確認
print("バンド幅 h^2: ", kde.covariance[0, 0])
# バンド幅 h^2:  0.9905809235665322

# Silvermanの因子が、定義通りの計算結果に等しいことを確認
print(kde.silverman_factor())
# 0.48438841363348917
print((kde.n * (kde.d + 2) / 4.)**(-1. / (kde.d + 4)))
# 0.4843884136334891

# データの普遍分散に Silvermanの因子の2乗を掛けてみると、先ほどのバンド幅の２乗に等しい
print((kde.silverman_factor()**2) * np.var(data, ddof=1))
# 0.9905809235665322

ちなみに、 statsmodelsのドキュメントを見てみましょう。
statsmodels.nonparametric.kde.KDEUnivariate.fit

The bandwidth to use. Choices are:
“scott” – 1.059 * A * nobs ** (-1/5.), where A is min(std(X),IQR/1.34)
“silverman” – .9 * A * nobs ** (-1/5.), where A is min(std(X),IQR/1.34)
“normal_reference” – C * A * nobs ** (-1/5.), where C is calculated from the kernel. Equivalent (up to 2 dp) to the “scott” bandwidth for gaussian kernels. See bandwidths.py
If a float is given, it is the bandwidth.

この様に、 “scott” のほうが、 $1.059*n^{(-0.2)}$ を掛けています。SciPyでは”silverman”の時に登場するものですね。
また、データの標準偏差よりも、IQR (これは 第3四分位数から第1四分位数を引いた値です)を1.34で割ったものが小さかったらそっちを使うと言うルールも実装されています。
SciPyのコードを見るとこれも入ってないです。

ちょっと色々疑わしいところはあるのですが、ぶっちゃけた話、ヒストグラムなどを使って可視化する時におまけにつけることがあるくらいのものなので、
あんまり気にしないで使って行こうかなと思っています。

ただ、statsmodelsの方が正しい値出してくれそうなのでこれもそのうちコードを載せて紹介しましょう。
次の記事ではガウスカーネル以外を使ったカーネル密度推定を先にやりたいので、Scikit-learnを使った実装を紹介する予定です。

前回の記事に続いてカーネル密度推定の話です。
参考: カーネル密度推定法の紹介

前の記事では数式に沿って実装するため、numpyで書きましたが、普段の利用では、SciPyにクラスが実装されているのでこれを使えます。

参考: scipy.stats.gaussian_kde
クラス名から分かる通り、カーネルはガウスカーネルがあらかじめ指定されており、他の関数は使えません。
実際のところ、ガウスカーネル以外が使えなくて困ることは滅多にないでしょう。

使い方ですが、 インスタンスを作るときに、データセットを渡すだけで完了します。
とてもありがたいことに、バンド幅の指定が要りません。
デフォルトでは Scottのルールと呼ばれる方法で決定されます。
bw_methodという引数で指定することで、もっと直接的な指定や、Silvermanのルールを使うことも可能です。
ScottのルールとSilvermanのルールがそれぞれ具体的にどの様にバンド幅を決めているのかは後日のブログで紹介するとして、
とりあえず、正規分布から生成したデータで使ってみましょう。

from scipy.stats import gaussian_kde
from scipy.stats import norm
import numpy as np
import matplotlib.pyplot as plt

# 50件のデータをサンプリング
data = norm(loc=2, scale=2).rvs(50)

# グラフ描写のため真の確率密度関数の値を取得しておく
xticks = np.linspace(-5, 9, 141)
pdf_values = norm(loc=2, scale=2).pdf(xticks)

# カーネル密度推定
kde = gaussian_kde(data)
# evaluate には pdfというエイリアスがあるのでそちらでも可
estimated_value = kde.evaluate(xticks)

# 可視化
fig = plt.figure(facecolor="w")
ax = fig.add_subplot(111)
ax.plot(xticks, pdf_values, label="真の分布")
ax.plot(xticks, estimated_value, label="推定値",  linestyle="--")
# 標本データ
for x in data:
    ax.vlines(x, ymin=0, ymax=0.01, alpha=0.5)
ax.set_ylim((0, 0.25))
ax.legend()
plt.show()

出力結果がこちらです。

とても手軽に良い推定が得られていますね。

ちなみに、推定してくれたバンド幅ですが、covariance という属性にその情報を持っています。

# バンド幅の確認 (hの2乗なので注意)
print(kde.covariance)

# [[0.79905207]]

多次元にも対応するため、分散共分散行列の形でデータを持てる様になっています。
今回は1次元データなので、ちょっと面倒ですが、無駄に2次元配列になってる中から要素を持ってくる必要があります。
varianceなので、この値は分散、前回の記事のバンド幅は標準偏差なので、平方根を取る必要があることも注意です。

せっかく推定してくれたバンド幅が取れるので、前回スクラッチで描いたコードの推定値と結果が一致することも見ておきましょう。
前回の記事のコピーだと、関数名が重複してしまって面倒なので適当に、my_kdeにリネームしました。

# ガウスカーネル
def gaussian_kernel(u):
    return(np.exp(-(u**2)/2)/np.sqrt(2*np.pi))


# ガウスカーネルを用いたカーネル密度推定
def my_kde(x, h, dataset):
    return sum([gaussian_kernel((x-d)/h) for d in dataset])/(len(dataset)*h)


# バンド幅の取得
h = np.sqrt(kde.covariance[0][0])
# 可視化
fig = plt.figure(facecolor="w")
ax = fig.add_subplot(111)
ax.plot(xticks, pdf_values, label="真の分布")
ax.plot(xticks, kde.evaluate(xticks), label="scipy", linestyle="--")
ax.plot(xticks, my_kde(xticks, h, data), label="スクラッチ", linestyle=":")
ax.legend()
plt.show()

SciPy で出した結果と、スクラッチで書いてるやつの結果が(一致してることの確認なので)重なってしまっていて、見にくいですが、
一応線の種類を変えてあるので、よく見るとしっかり同じところに線が引かれているのが分かります。

今回からしばらくカーネル密度推定法について書いていきます。

参考文献ですが「パターン認識と機械学習 上 (PRML 上)」の 119ページ、 2.5.1 カーネル密度推定法　あたりがわかりやすいのではないでしょうか。
ただし、この本だと最後にパラメーター$h$の決め方が重要だ、ってところまで書いて、$h$の決め方を紹介せずに次の話題に移っているので、
これからの数回の記事の中でそこまで説明できたらいいなと思っています。

さて、カーネル密度推定法を大雑把に言うと、
「確率密度が全然わからない確率分布から標本が得られたときに、その標本からできるだけいい感じに元の確率密度を推定する方法」です。
元の確率密度が正規分布だとかポアソン分布だとかだけでもわかっていたら、その期待値や分散などのパラメーターだけ推定すればいい(これをパラメトリックという)のですが、
それさえもわからないときに使われる手法であり、ノンパラメトリックなアプローチと呼ばれます。

できるだけ一般的な形でカーネル密度推定法の方法を書くと次のような形になります。

ある$D$次元のユークリッド空間の未知の確率密度$p(\mathbf{x})$から、観測値の集合$\{\mathbf{x}_n\}\ \ (n=1,2, \cdots, N)$が得られたとします。

そして、次の性質を満たすカーネル関数$k(\mathbf{u})$を選びます。
$$
\begin{align}
k(\mathbf{u}) &\geq 0,\\
\int k(\mathbf{u}) d\mathbf{u} &= 1
\end{align}
$$
このとき、$p(\mathbf{x})$を次で推定するのがカーネル密度推定法です。($h$はパラメーター。)
$$
p(\mathbf{x}) = \frac{1}{N}\sum_{n=1}^{N}\frac{1}{h^D}k(\frac{\mathbf{x}-\mathbf{x}_n}{h}).
$$

カーネル密度関数は要は確率密度間数ですね。理論上は何でもいいそうですが、現実的には原点周りで値が大きく、そこから外れると0に近く（もしくは等しく）なる関数が好都合です。
正規分布の確率密度関数(ガウスカーネル)が頻繁に用いられます。

上記の式が導出されるまでの話はPRMLに書いてあります。
その部分を省略してしまったので代わりに上の定義通り計算して、未知の確率分布が近似できている様子を確認しておきましょう。
できるだけライブラリ使わないでやります。(ただ、正解データの分布作成でscipyは使っていますし、numpyやmatplotlibは例外とします。)
カーネル関数にはガウスカーネルと、単位超立方体内で1、その外で0を取る関数(パルツェン窓というらしい)の2種類でやってみます。

元となる分布は正規分布２個を足した、２つ山のある分布で実験します。($D$は1です)

最初に正解になる分布関数の定義と、そこからのデータサンプリングを済ませておきます。

from scipy.stats import norm
import numpy as np
import matplotlib.pyplot as plt


# 推定したい分布の真の確率密度関数
def p_pdf(x):
    return (norm.pdf(x, loc=-4, scale=2)+norm.pdf(x, loc=2, scale=1))/2


# その分布からのサンプリング
def p_rvs(size=1):
    result = []
    for _ in range(size):
        if norm.rvs() < 0:
            result.append(norm.rvs(loc=-4, scale=2))
        else:
            result.append(norm.rvs(loc=2, scale=1))
    return np.array(result)


# グラフ描写のため真の確率密度関数の値を取得しておく
xticks = np.linspace(-10, 5, 151)
pdf_values = p_pdf(xticks)

# 50件のデータをサンプリング
X = p_rvs(50)

そして、それぞれのカーネル関数と、それを使ったカーネル密度関数を定義しておきます。

# パルツェン窓
def parzen_window(u):
    if np.abs(u) <= 0.5:
        return 1
    else:
        return 0


# ガウスカーネル
def gaussian_kernel(u):
    return(np.exp(-(u**2)/2)/np.sqrt(2*np.pi))


# パルツェン窓を用いたカーネル密度推定
def parzen_kde(x, h, dataset):
    return sum([parzen_window((x-d)/h) for d in dataset])/(len(dataset)*h)


# ガウスカーネルを用いたカーネル密度推定
def gaussian_kde(x, h, dataset):
    return sum([gaussian_kernel((x-d)/h) for d in dataset])/(len(dataset)*h)

あとは可視化です。 今回は パラメーター$h$は$0.5, 1, 2$の３つ試しました。

fig = plt.figure(figsize=(12, 12), facecolor="w")

# 3種類のhで実験
for i, h in enumerate([0.5, 1, 2]):
    parzen_kde_values = [parzen_kde(x, h, X) for x in xticks]
    gaussian_kde_values = [gaussian_kde(x, h, X) for x in xticks]

    ax = fig.add_subplot(3, 2, i*2+1, facecolor="w")
    ax.set_title(f"パルツェン窓, h={h}")
    ax.plot(xticks, pdf_values, label="真の分布")
    ax.plot(xticks, parzen_kde_values, label="パルツェン窓")
    # 標本データ
    for x in X:
        ax.vlines(x, ymin=0, ymax=0.01)
    ax.set_ylim((0, 0.37))
    ax.legend()

    ax = fig.add_subplot(3, 2, i*2+2, facecolor="w")
    ax.set_title(f"ガウスカーネル, h={h}")
    ax.plot(xticks, pdf_values, label="真の分布")
    ax.plot(xticks, gaussian_kde_values, label="ガウスカーネル")
    # 標本データ
    for x in X:
        ax.vlines(x, ymin=0, ymax=0.01)
    ax.set_ylim((0, 0.37))
    ax.legend()

plt.show()

出力結果がこれです。

青い線が真の分布、黒くて短い縦線たちが推定に使った標本データ、
オレンジの線が推定した分布です。

「元の確率密度関数の形が全くわからない」という前提から出発している割にはうまく推定できているのもあると言えるのではないでしょうか。
また、カーネル密度関数の選び方、$h$の決め方がそれぞれ重要であることも分かりますね。

やってることといえば、得られた標本の付近で大きめの値を取る関数を全ての標本について計算してその平均を取っているだけなのですが、シンプルな割に強力です。