シンプソンのパラドックス

先日、データを分析している中でシンプソンのパラドックスが発生しているのを見かけました。
興味深い現象なので、紹介したいと思います。
ただし、業務的な情報は書けないので記事中の用語も設定もデータも全部架空の物です。

2種類のアプリがあったとします。それぞれ旧アプリと、新アプリとします。
そしてそれらのアプリを使っているユーザーがとある属性によってグループA,グループBに分かれていたとします。

ユーザー数の内訳が次のようになっていたとします。(単位:人)

旧アプリ 新アプリ
グループA 40000 1000
グループB 60000 9000

これらのユーザーのコンバージョン数が次の通りだったとします。

旧アプリ 新アプリ
グループA 3200 100
グループB 1800 360

コンバージョン率を見ると次のようになりますね。

旧アプリ 新アプリ
グループA 8% 10%
グループB 3% 4%

どちらのグループのユーザーに対しても、新アプリの方がコンバージョン率が高いことがわかりました。

しかしここで、グループごとに分けて集計することをやめて、新アプリと旧アプリを単純に比較してみます。

旧アプリ 新アプリ
ユーザー数 100000 10000
コンバージョン数 5000 460
コンバージョン率 5% 4.6%

なんと、新アプリより旧アプリの方がコンバージョン率が高いことになりました。

このように、
集団全体を複数の集団に分けてそれぞれの集団で同じ仮説(今回の例では新アプリの方がコンバージョン率が高い)が成り立っても、
集団全体に対してはそれが成り立たないことがあることをシンプソンのパラドックスと呼びます。

これは$\frac{a}{A}>\frac{b}{B}$ かつ $\frac{c}{C}>\frac{d}{D}$ が成り立ったとしても
$$\frac{a+c}{A+C}>\frac{b+d}{B+D}$$
が成り立つわけじゃないと言う単純な数学的な事実から発生する物です。

今回の例で言えば、新アプリの方がどちらのユーザー群に対しても良い効果をもたらしているので良さそうなのに、
全体の集計だけで旧アプリの方が良いと結論づけてしまうと誤った分析をしてしまうことになります。
注意する必要がありますね。

SciPyで多次元のカーネル密度推定

以前も紹介した、SciPyのカーネル密度推定のメソッド、gaussian_kdeの話です。
参考: SciPyによるカーネル密度推定

最近、多次元(と言っても2次元のデータですが)に対して、カーネル密度推定を行いたいことがあり、
どうせ1次元の場合と同じように使えるのだろうと適当に書いたら思うような動きにならず苦戦しました。

何をやろうとしたかというと、
[[x0, y0], [x1, y1], [x2, y2], …, [xn, yn]]
みたいなデータをそのままgaussian_kdeに渡してしまっていました。

ドキュメント: scipy.stats.gaussian_kde
をよく読むと、
Datapoints to estimate from. In case of univariate data this is a 1-D array, otherwise a 2-D array with shape (# of dims, # of data).
と書いてあります。

僕が渡そうとしたデータは shapeが[データ件数, 2(=次元数)]になっていたのですが、実際は
[2(=次元数), データ件数]の型で渡さないといけなかったわけです。

丁寧なことに、Examplesに取り上げられている例も1Dではなく2Dの例で、np.vstackとか使って書かれているので、以前の記事を書いた時にもっとしっかり読んでおけばよかったです。

使い方がわかったので、2次元のデータに対してやってみます。
サンプルデータは今回はscikit-learnのmake_moonsを使いました。


from sklearn.datasets import make_moons
import matplotlib.pyplot as plt

# データ生成。 今回はラベルは不要なのでデータだけ取得する。
data, _ = make_moons(
    n_samples=200,
    noise=0.1,
    random_state=0
)

fig = plt.figure(facecolor="w")
ax = fig.add_subplot(111, title="サンプルデータ")
ax.scatter(data[:, 0], data[:, 1])
plt.show()

サンプルデータはこんな感じです。

生成したデータを、gaussian_kdeにそのまま渡すとうまくいきません。
出来上がるモデルは、200次元のデータ2個を学習したものになっているからです。


from scipy.stats import gaussian_kde
# これはエラーは出ないが誤り
kde = gaussian_kde(data)

# 密度推定した結果を得ようとするとエラーになる。
print(kde.pdf([0, 0]))
# ValueError: points have dimension 1, dataset has dimension 200

きちんと、2次元のデータ200個を学習させるには、データを転置させて渡します。


from scipy.stats import gaussian_kde
# (データを、(次元数, データ件数) の型で渡す)
kde = gaussian_kde(X.T)

kde.evaluate を使って、推定した結果の値を得る時も、(次元数, データ件数)の形で、データを渡す必要があります。
(一点のみなら長さが次元数分の配列を渡せば良いです。)


# 点(1, 1)での値
print(kde.evaluate([1, 1]))
# [0.03921808]
# 4点(-0.5, 0), (0, 1), (0.5, 1), (1, 0) での値
print(kde.evaluate([[-0.5, 0, 0.5, 1], [0, 1, 1, 0]]))
# [0.07139624 0.2690079  0.2134083  0.16500181]

等高線を引いて図示する場合は次のように行います。(公式ドキュメントでは等高線ではなく、imshowで可視化していますね。)
こちらの記事も参考にしてください。
参考: matplotlibで等高線


# 等高線を引く領域のx座標とy座標のリストを用意する
x = np.linspace(-1.5, 2.5, 41)
y = np.linspace(-0.8, 1.3, 22)
# メッシュに変換
xx, yy = np.meshgrid(x, y)
# kdeが受け取れる形に整形
meshdata = np.vstack([xx.ravel(), yy.ravel()])
# 高さのデータ計算
z = kde.evaluate(meshdata)

# 可視化
fig = plt.figure(facecolor="w")
ax = fig.add_subplot(111, title="カーネル密度推定")
ax.scatter(data[:, 0], data[:, 1], c="b")
ax.contourf(xx, yy, z.reshape(len(y), len(x)), cmap="Blues", alpha=0.5)
plt.show()

結果がこちらです。

最後に、xもしくはy座標を固定して、断面をみる方法を紹介しておきます。
これはシンプルに、固定したい方は定数値でもう一方のデータと同じ長さの配列を作って、
固定しない方のデータを動かしてプロットするだけです。
x=-1.0, 0.5, 2.0 の3つの直線で切ってみた断面を可視化すると次のようなコードになります。


fig = plt.figure(facecolor="w")
ax = fig.add_subplot(111)
for x_ in [-1.0, 0.5, 2.0]:
    ax.plot(y, kde.pdf([[x_]*len(y), y]), label=f"x={x_}")
ax.set_ylabel("z")
ax.set_xlabel("y")
ax.legend()
plt.show()

結果がこちら。

だいぶ使い方の感覚が掴めてきました。

Scipyで対数正規分布

対数正規分布の話の続きです。
参考: 前回の記事

Pythonで対数正規分布を扱う場合、(もちろんスクラッチで書いてもいいのですが普通は)Scipyに実装されている、
scipy.stats.lognormを使います。
これに少し癖があり、最初は少してこずりました。

まず、対数正規分布の確率密度関数はパラメーターを二つ持ちます。
次の式の$\mu$と$\sigma$です。
$$
f(x) =\frac{1}{\sqrt{2\pi}\sigma x}\exp\left(-\frac{(\log{x}-\mu)^2}{2\sigma^2}\right).
$$

それに対して、scipy.stats.lognorm の確率密度関数(pdf)は、
s, loc, scale という3つのパラメーターを持ちます。 ややこしいですね。

正規分布(scipy.stats.norm)の場合は、 loc が $\mu$に対応し、scaleが$\sigma$に対応するので、とてもわかりやすいのですが、lognormはそうはなっていなっておらず、わかりにくくなっています。
(とはいえ、ドキュメントには正確に書いてあります。)

Scipyの対数正規分布においては、確率密度関数は
$$
f(x, s) = \frac{1}{sx\sqrt{2\pi}}\exp\left(-\frac{\log^2{x}}{2s^2}\right)
$$
と定義されています。
そして、$y=(x-loc)/scale$と変換したとき、lognorm.pdf(x, s, loc, scale) は、 lognorm.pdf(y, s) / scale と等しくなります。(とドキュメントに書いてあります。)
わかりにくいですね。

$\mu$と$\sigma$とはどのように対応しているのか気になるのですが、上の2式を見比べてわかるとおり、(そしてドキュメントにも書いてある通り、)
まず、引数のsが、パラメーターの$\sigma$に対応します。 (scaleじゃないんだ。)
そして、引数のscaleは$e^{\mu}$になります。(正規分布の流れで、locと$\mu$が関係してると予想していたのでこれも意外です。)
逆にいうと、$\mu=\log{scale}$です。

locはまだ登場していませんが、一旦ここまでの内容をプログラムで確認しておきます。
lognorm(s=0.5, scale=100) とすると、 $\sigma=0.5$で、$\mu=\log{100}\fallingdotseq 4.605$の対数正規分布になります。
乱数をたくさん取って、その対数が、期待値$4.605…$、標準偏差$0.5$の正規分布に従っていそうかどうかみてみます。


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

data = lognorm(s=0.5, scale=100).rvs(size=10000)
log_data = np.log(data)
print("対数の平均:", log_data.mean())
# 対数の平均: 4.607122120944554
print("対数の標準偏差:", log_data.std())
# 対数の標準偏差: 0.49570170767616645

fig = plt.figure(figsize=(12, 6), facecolor="w")
ax = fig.add_subplot(1, 2, 1, title="元のデータ")
ax.hist(data, bins=100, density=True)
ax = fig.add_subplot(1, 2, 2, title="対数")
ax.hist(log_data, bins=100, density=True)
# 期待値のところに縦線引く
ax.vlines(np.log(100), 0, 0.8)

plt.show()

出力された図がこちらです。

想定通り動いてくれていますね。

残る loc ですが、 これは分布の値を左右に平行移動させるものになります。
先出の$y=(x-loc)/scale$ を $x = y*scale + loc$ と変形するとわかりやすいかもしれません、

わかりやすくするために、locに極端な値($\pm 1000$と$0$)を設定して、乱数をとってみました。
(sとscaleは先ほどと同じです。)


locs = [-1000, 0, 1000]

fig = plt.figure(figsize=(18, 6), facecolor="w")
for i, loc in enumerate(locs, start=1):
    data = lognorm(s=0.5, scale=100, loc=loc).rvs(size=1000)
    ax = fig.add_subplot(1, 3, i, title=f"loc={loc}")
    ax.hist(data, bins=100, density=True)

plt.show()

出力がこちらです。

分布の左端が それぞれ、locで指定した、 -1000, 0, 1000 付近になっているのがみて取れると思います。

対数正規分布

最近とあるデータのモデリングのために対数正規分布を使うことがありました。
そのとき使ったScipyのlognormには少し癖があったのでそれを紹介しようと思ったのですが、その前に対数正規分布そのものについて紹介させてもらいます。

対数正規分布とは一言で言ってしまうと、「それに従う確率変数の対数を取ったら正規分布に従うような確率分布」です。
正規分布の対数ではなく、対数を取ったら正規分布なので、注意が必要です。

もう少し正確に書きます。
確率変数$Y$が、正規分布$N(\mu, \sigma^2)$に従うとします。
このとき、$X=e^{Y}$は対数正規分布に従います。

定義に従って確率密度関数$f(x)$を導出しておきましょう。
まず、正規分布 $N(\mu, \sigma^2)$ の確率密度間数$g(y)$は、
$$
g(y) = \frac{1}{\sqrt{2\pi}\sigma}\exp\left(-\frac{(y-\mu)^2}{2\sigma^2}\right)
$$
です。
あとは$y=\log{x}$と$\frac{dy}{dx}=\frac{1}{x}$ に気をつけて、次のように計算できます。

$$
\begin{align}
f(x) &= g(y)\frac{dy}{dx}\\
&=g(\log{x})\frac{1}{x}\\
\therefore f(x) &=\frac{1}{\sqrt{2\pi}\sigma x}\exp\left(-\frac{(\log{x}-\mu)^2}{2\sigma^2}\right).
\end{align}
$$

ここで、$x$の値域は、 $0 <x < \infty$ です。
$\mu$ と $\sigma$ はパラメーターになります。
この導出の過程から明らかですが、それぞれ、対数正規分布に従う確率変数の「対数の」期待値と標準偏差です。
この二つがそのまま対数正規分布の期待値や標準偏差になるわけではないので気をつける必要があります。

対数ではない、確率変数$X$自体の期待値と分散は次の値になります。
参考: 対数正規分布 – Wikipedia
$$
\begin{align}
E[X] &= e^{\mu+\sigma^2/2}\\
V[X] &= e^{2\mu+\sigma^2}(e^{\sigma^2}-1)
\end{align}
$$

バートレット検定

前回の記事の中で、SciPyにはF検定の関数が実装されていないという話を書きました。
参考: PythonでF検定を実装する

しかし、その代わりにSciPyには分散が等しいことを検定する別の方法が実装されています。
その一つが、バートレット検定です。
関数はこれです。
scipy.stats.bartlett

複数(2個だけでなく3個以上も可能)のサンプルを渡してあげると、統計量とp値を計算してくれます。

今回の記事では、使い方の前にバートレット検定の統計量を紹介します。
参考にしたのは、東京大学出版会の自然科学の統計学 (青い本) の 94ページです。

サンプルの数を$a$とし、それぞれのサンプルのサイズを$n_i \ \ (j=1,\dots, a)$とします。
また、サンプルサイズの合計を$n=\sum_{j}n_j$と置いておきます。
各サンプルの不偏分散を
$$
V_i = \sum_{j} \frac{(y_{ij}-\bar{y_i})^2}{n_i-1}
$$
として、さらにこれらを併合したものを
$$
V_e = \sum_i \frac{(n_i-1)V_i}{n-a} = \sum_i \sum_j \frac{(y_{ij}-\bar{y_i})^2}{n-a}
$$
とします。
さらに、
$$
B = (n-a)\log{V_e} -\sum_i (n_i-1)\log{V_i}
$$
と置いて、
$$
B’ = \frac{B}{
1+\frac{1}{3(a-1)}\{\sum_i \frac{1}{n_i-1} – \frac{1}{n-a}\}
}
$$
と補正すると、この$B’$は自由度$a-1$の$\chi^2$分布に近似的にしたがいます。
それを利用して、検定を行うのがバートレット検定です。

数式に沿ってNumPyで実装してみたのが次の関数です。ランダム生成した3サンプルで試してみました。


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

# スクラッチで実装したバーレット検定量の算出関数
def bartlett_value(*args):
    # サンプル数
    a = len(args)
    # 各サンプルの サンプルサイズと不偏分散
    ni = np.zeros(a)
    Vi = np.zeros(a)
    for j in range(a):
        ni[j] = len(args[j])
        Vi[j] = np.var(args[j], ddof=1)

    # サンプルサイズの合計
    N = np.sum(ni)
    # Veの計算
    Ve = np.sum((ni-1)*Vi) / (N-a)

    # 統計量の計算
    B = (N-a) * np.log(Ve) - np.sum((ni-1)*np.log(Vi))
    C = 1+1/(3*(a-1))*(np.sum(1/(ni-1))-1/(N-a))

    return B/C


# サンプルを3種類生成する
data_a = norm(loc=3, scale=5).rvs(20)
data_b = norm(loc=7, scale=8).rvs(30)
data_c = norm(loc=-5, scale=5).rvs(15)

# 統計量を算出
b_value = bartlett_value(data_a, data_b, data_c)

# p値を算出
p_value = chi2(3-1).sf(b_value)

print("B値:", b_value)
# B値: 8.26319721031961
print("p値:", p_value)
# p値: 0.016057189188779606

$\alpha=0.05$とすると帰無仮説が棄却され、分散が等しくないことが検出されていますね。

もちろん、SciPyに関数は用意されているので、普段はこの様にスクラッチで実装する必要はなく、
呼び出すだけで良いです。結果が一致することを見ておきます。


from scipy.stats import bartlett


b_value, p_value = bartlett(data_a, data_b, data_c)
print("B値:", b_value)
# B値: 8.26319721031961
print("p値:", p_value)
# p値: 0.016057189188779606

バッチリ一致しました。

statsmodels で カーネル密度推定

今回は 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でカーネル密度推定

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

先ほどの 6種類のカーネル + 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のガウシアンカーネル密度推定に実装されているバンド幅の算出方法について

僕の環境の 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のコードもかなり参考になりました。

さて、推定したバンド幅はその2乗の値が covariance という属性に格納されています。(1次元データの場合。)
これは、 学習したデータの不変分散に、 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乗を掛けてみると、先ほどのバンド幅の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乗を掛けてみると、先ほどのバンド幅の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を使った実装を紹介する予定です。

SciPyによるカーネル密度推定

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

前の記事では数式に沿って実装するため、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種類でやってみます。

元となる分布は正規分布2個を足した、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$の3つ試しました。


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$の決め方がそれぞれ重要であることも分かりますね。

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