ガンマ分布がポアソン分布の共役事前分布である事の証明

前回の記事に続いて今回もガンマ分布の話です。そもそもこのブログの一連の記事でガンマ分布が登場した理由はこれがポアソン分布の共役事前分布だからですが、有名な性質の割にこれの証明を自分がちゃんとやったことがないような気がしたのでやってみました。

共役事前分布とは

最初に軽く共役事前分布自体の説明を書いておきます。これはベイズ統計における概念で、特定の統計モデルにおいて、事後分布の形が事前分布と同じ分布族に属するような事前分布のことを指します。

例えば、ベータ分布は二項分布の共役事前分布なので、二項分布のパラメーター$p$の事前分布をベータ分布にしておくと、事後分布もベータ分布になるので計算が楽です。

早速本題のガンマ分布がポアソン分布の共役事前分布であることを証明してみましょう。

ポアソン分布の尤度関数

まず、ポアソン分布の確率密度関数は期待値$\lambda$をパラメーターとして以下になります。

$$
Po(k|\lambda) = \frac{e^{-\lambda}\lambda^k}{k!}.
$$

ここで観測値$\{x_1, x_2, _dots, x_n\}$が与えられると、その尤度関数は次になります。

$$
L(\lambda) = \prod_{i=1}^{n} Po(x_i|\lambda) = \prod_{i=1}^{n} \frac{e^{-\lambda}\lambda^{x_i}}{x_{i}!} = e^{-n\lambda}\lambda^{\sum_{i=1}^n x_i} \prod_{i=1}^{n} \frac{1}{x_{i}!}.
$$

事前分布のガンマ分布

次に事前分布を見ておきます。以前の記事で$k$, $\lambda$をパラメーターとして使いましたがこの記事では先に別の意味でこれらの文字を使ってしまっていますので$\alpha$, $\beta$をパラメーターとして使うと、$\lambda$の事前分布の確率密度関数は次で表されます。

$$
Ga(\lambda| \alpha, \beta) = \frac{\beta^\alpha}{\Gamma(\alpha)}\lambda^{\alpha-1}e^{-\beta\lambda}.
$$

事後分布の導出

さて、ここから事後分布を計算していきましょう。分布の形だけ着目すればいいので、事後分布の$\lambda$に関係ない部分は無視していきます。

$$\begin{align}
p(\lambda| x_1, x_2, _dots, x_n) & \propto L(\lambda) \times Ga(\lambda| \alpha, \beta)\\
& \propto \left(e^{-n\lambda}\lambda^{\sum_{i=1}^n x_i} \right) \times \left( \lambda^{\alpha-1}e^{-\beta\lambda} \right)\\
& \propto \lambda^{\alpha -1 +\sum_{i=1}^{n} x_i} e^{-\lambda(\beta+n)}.
\end{align}
$$

さて、この導出された事後分布(の確率密度関数に比例する式)の形式をよく見るとガンマ分布の形をしていることがわかります。具体的には、次のようになります。

$$
p(\lambda| x_1, x_2, _dots, x_n) = Gamma(\lambda | \alpha +\sum_{i=1}^{n} x_i, \beta+n)
$$

以上で、事後分布もガンマ分布であり、ガンマ分布がポアソン分布の協約事前分布であることが証明されました。

ガンマ分布について

前々回の記事の中で例として挙げたモデルの中でガンマ分布を使ったのですが、そういえばガンマ分布を久しぶりに使ったなぁと思ったのでこの機会にその定義や性質等をまとめておこおうと思います。
参考: pyMCで変化点検出 (事前分布としてガンマ分布を採用しました)

まず、ガンマ分布の名前の由来になっているのがガンマ関数なので先にその定義も書いておきます。ガンマ関数は実部が正の複素数$z$に対して次の積分で定義される関数です。ちなみに解析接続して、一般の複素数に対しても定義されます。

$$
\Gamma(z) = \int_{0}^{\infty} t^{z-1}e^{-t}dt.
$$

このガンマ関数は、$$\Gamma(n+1) = n!$$という極めて重要な関係式が成り立ち、階乗の一般化とみなせます。

本題のガンマ関数ですが、これは形状母数$k>0$と尺度母数$\theta>0$を用いて確率密度関数が次の式で定義される関数です。

$$
f(x) = \frac{1}{\Gamma(k)\theta^k}x^{k-1}e^{-x/\theta} \qquad (x > 0) .
$$

$\theta$の代わりに、$\lambda = \frac{1}{\theta}$を用いて、次のように表現されることもあります。僕は個人的にはこっちの方が好きです。

$$
f(x) = \frac{\lambda^k}{\Gamma(k)}x^{k-1}e^{-\lambda x} \qquad (x > 0) .
$$

期待値と分散は次の式で表されます。後で紹介しますが、指数分布との関係を考慮するととても自然な結果です。

$$
\begin{align}
E(X)&=k\theta&=\frac{k}{\lambda}\\
V(X)&=k\theta^2&=\frac{k}{\lambda^2}\\
\end{align}
$$

他の確率分布との関係

ガンマ分布はさまざまな分布との関係が深い確率分布です。特に基本的なものを紹介します。

まず、ガンマ分布の$k$が整数の場合に特別な名前がついていて、これをアーラン分布と言います。ガンマ関数を階乗として書けますので、次の形になりますね。
$$f(x) = \frac{1}{(k-1)!\theta^k}x^{k-1}e^{-x/\theta} \qquad (x > 0) .$$

個人的に重要だと考えているのは、互いに独立で同一の指数分布に従う$k$個の確率変数の和は形状母数$k$のガンマ分布に従うことです。個数$k$は明らかに整数なのでこれは正確にはアーラン分布の性質なのですが、ガンマ分布の特徴と考えても良いでしょう。

逆にガンマ分布において$k=1$とすると、ただの指数分布になります。

指数分布についてはこのブログでも過去に紹介しています。
参考: 指数分布について

その記事で紹介した指数分布の期待値と分散を考慮すると、ガンマ分布の期待値と分散がちょうど$k$倍になるのは納得ですね。

もう一つ、1以上の整数$n$に対して、$k=\frac{n}{2}$かつ$\theta=2$とすると、ガンマ分布は自由度$n$のカイ二乗分布になります。
カイ二乗分布もこのブログで紹介しているので確率密度関数の式の形を見比べてみてください。
参考: カイ二乗分布について

そして最後に、ガンマ分布はポアソン分布とも深い関係にあります。ベイズ統計においてはガンマ分布がポアソン分布の共役事前分布になることが特に重要でしょう。

前出の変化点検出のモデルでも、この性質があるのでガンマ分布を事前分布に採用しました。これについては別記事でちゃんと証明をつけたいと思います。

scipy.statsで確率分布を自作する

SciPyのstatsモジュールには非常に多くの確率分布が定義されています。
参考: Statistical functions (scipy.stats) — SciPy v1.12.0 Manual

ほとんどの用途はこれらを利用すると事足りるのですが、自分で定義した確率分布を扱うこともあり、scipyのstatsに用意されているような確率分布と同じようにcdfとかのメソッドを使いたいなと思うことがあります。

そのような場合、SciPyではscipy.stats.rv_continuousを継承してpdfメソッドを定義することで確率分布を定義でき、そうすれば(計算量や時間の問題などありますが)cdf等々のメソッドが使えるようになります。

参考: scipy.stats.rv_continuous — SciPy v1.12.0 Manual

実はSciPyの元々用意されている確率分布たちもこのrv_continuousを継承して作られているので、それらのソースコードを読むと使い方の参考になります。既存の確率分布等はpdfだけでなく、cdf、sf、ppf、期待値や分散などの統計量などが明示的に実装されていて計算効率よく使えるようになっています。

しかし先ほども述べた通り、pdf以外は必須ではなく実装されていなければpdfから数値的に計算してくれます。(ただ、予想外のところで計算が終わらなかったり精度が不十分だったりといった問題が起きるので可能な限り実装することをお勧めします。)

確率分布クラスの自作方法

さて、それでは実際にやってみましょう。既存に存在しない例が良いと思うので、正規分布を2個組み合わせた、山が2個ある分布を作ってみます。

クラスを継承して _pdf (pdfではない) メソッドをオーバーライドする形で実装します。

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


# 確率密度関数の定義
def my_pdf(x):
    curve1 = np.exp(-(x+3)**2 / 2) / np.sqrt(2 * np.pi)
    curve2 = np.exp(-(x-3)**2 / 2) / np.sqrt(2 * np.pi)

    return (curve1 + curve2)/2


# rv_continuousクラスのサブクラス化
class MyDistribution(rv_continuous):
    def _pdf(self, x):
        return my_pdf(x)


# 分布のインスタンス化
my_dist = MyDistribution(name='my_dist')


# 可視化
x = np.linspace(-6, 6, 121)
fig = plt.figure(facecolor="w", figsize=(8, 8))
ax = fig.add_subplot(2, 2, 1, title="pdf")
ax.plot(x, my_dist.pdf(x))

ax = fig.add_subplot(2, 2, 2, title="cdf")
ax.plot(x, my_dist.cdf(x))

ax = fig.add_subplot(2, 2, 3, title="sf")
ax.plot(x, my_dist.sf(x))

ax = fig.add_subplot(2, 2, 4, title="logpdf")
ax.plot(x, my_dist.logpdf(x))

plt.show()

結果がこちらです。

実装したのはpdfだけでしたが、cdfやsf、logpdfなども動いていますね。
繰り返しになりますが、これらは数値計算されたものなので十分注意して扱ってください。特に精度と計算量はもちろんですが、指数関数等もからむので、極端な値を入れたら計算途中でオーバーフローが起きたりもします。重要なものは_pdfと同様に_cdfなどとして明示的に定義した方が良いでしょう。

たとえば、先ほどの分布は分散の計算に少し時間がかかります。

%%time
my_dist.var()
"""
CPU times: user 8.02 s, sys: 146 ms, total: 8.16 s
Wall time: 8.06 s
10.000000000076088
"""

確率密度関数を用意するときの注意点

確率密度関数を自分で作りましたが、この内容に対してのバリデーションなどは用意されていないので自分で責任を持って管理する必要があります。

たとえば、変数xの定義域で正であること、定義域全体で積分したら1になることなどは事前に確認しておく必要があります。(確率密度関数の定義を満たしてないととcdfなど他のメソッドが変な動きをします)

確率分布の台が有限の場合

もう一点、先ほどの2山の分布はxが$-\infty$から$\infty$の値を取る分布でしたが、確率分布の中にはそうではないものもあります。正の範囲でしか定義されない対数正規分布や指数分布、有限区間でしか定義されないベータ分布などですね。

これらについては、「分布のインスタンス化」を行うタイミングで台の下限上限をそれぞれa, b という引数で行います。省略した方は$-\infty$もしくは$\infty$となります。

確率密度関数が上に凸な放物線で定義された分布でやるとこんな感じです。

def my_pdf2(x):
    return 6*x*(1-x)


class MyDistribution2(rv_continuous):
    def _pdf(self, x):
        return my_pdf2(x)


# 分布のインスタンス化
my_dist2 = MyDistribution2(a=0, b=1, name='my_dist2')

以上が、SciPyで連続確率分布インスタンスを自作する方法でした。

ジニ不純度について

ここ数週間にわたってエントロピー関連の話が続きましたので、ついでみたいになるのですが決定木の評価関数として同じようによく使われるジニ不純度についても紹介しておきます。
ジニ係数(はじパタなど)や、ジニ指数(カステラ本など)といった呼び名もあるのですが、経済学で所得格差の指標に用いられるジニ係数と紛らわしくなるので、この記事ではジニ不純度と呼びます。

ジニ不純度の定義

ジニ不純度は、決定木やランダムフォレストなどの機械学習のアルゴリズムにおいて、ノードの純度を計測するための指標の一つとして用いられます。要するに、そのノードに分類される観測値のクラスの内訳が全部同じだとなると不純度は小さく、同じノードに複数のクラスの観測値が混在してたら不純度は高くなります。エントロピーと似た挙動ですね。

数式としては次のように定義されます。$n$クラスに分類する問題の場合、$i$番目のクラスのサンプル比率を$p_i$とすると、ジニ不純度は次のように定義されます。

$$I = 1 – \sum_{i=1}^n p_i^2.$$

例えば、3クラスの分類問題を考えると、あるノードの観測値が全部同じクラスだったとします。すると、$p_1 = 1, p_2=0, p_3 = 0$みたいになるのですが、この時$I=0$と、不純度が0になります。また、3クラスが均等に混ざっていた場合、$p_1 = p_2 = p_3 = 2/3$となり、時に不純度は$I = 2/3$となります。$n$クラスの分類問題の場合、この$1-1/n$が最大値です。

もう7年近く前になりますが、初めてこの定義式を見た時は、1から確率の平方和を引くことに対してちょっと違和感を感じました。ただ、これは計算を効率的にやるためにこの形にしているだけで、実際は次の形で定義される値と考えておくと納得しやすいです。

$$I = \sum_{i=1}^n p_i (1-p_i).$$

これ展開すると、$\sum_{i=1}^n p_i – p_i^2$となり、$\sum_{i=1}^n p_i =1$ですから、先に出した定義式と一致することがわかります。

ジニ不純度の解釈

それでは、$p_i(1-p_i)$とは何か、という話なのですが解釈が二つあります。

まず一つ目。ある観測値がクラスiであれば1、クラスi出なければ0というベルヌーイ試行を考えるとその分散がこの$p_i(1-p_i)$になります。これを各クラス分足し合わせたのがジニ不純度ですね。

クラスiである確率が0の場合と1の2種類の場合がまじりっけ無しで不純度0ということでとてもわかりやすいです。

もう一つは、そのノードにおける誤分類率です。あるノードにおける観測値を、それぞれ確率$p_i$で、クラス$i$である、とジャッジすると、それが本当はクラス$i$ではない確率が$(1-p_i)$ありますから、誤分類率は$\sum_{i=1}^n p_i (1-p_i)$となります。

あるノードの観測値のクラスが全部一緒なら誤分類は発生しませんし、$n$クラスが均等に含まれていたら誤分類率は最大になりますね。不純度の定義としてわかりやすいです。

エントロピーと比較した場合のメリット

決定木の分割の評価として使う場合、エントロピーとジニ不純度は少しだけ異なる振る舞いをし、一長一短あるのですが、ジニ不純度には計算効率の面で優位性があるとも聞きます。というのも、$p\log{p}$より$p^2$のほうが計算が簡単なんですよね。scikit-learnがginiの方をデフォルトにしているのはこれが理由じゃないかなと思ってます。(個人的な予想ですが。)

ただ、実際に実験してみると計算速度はそこまで大差はないです(np.log2の実装は優秀)。なんなら自分で書いて試したらエントロピーの方が速かったりもします。とはいえ、$p\log{p}$の計算はちゃんとやるなら、$p=0$ではないことをチェックして場合分けを入れたりしないといけないのでそこまできちんと作ったらジニ不純度の方が速くなりそうです。

Pythonでの実装

数式面の説明が続きましたが、最後に実装の話です。実は主要ライブラリではジニ不純度を計算する便利な関数等は提供されていません。まぁ、簡単なので自分で書けば良いのですが。

scikit-learnのソースコードを見ると、ここで実装されて内部で使われていますね。

確率$p$のリストが得られたら、全部二乗して足して1から引けば良いです。
$1/6+1/3+1/2=1$なので、確率のリストは$[1/6, 1/3, 1/2]$としてやってみます。

import numpy as np


p_list = np.array([1/6, 1/3, 1/2])
print(1-(p_list**2).sum())
# 0.6111111111111112

まぁ、確かにこれならわざわざライブラリ用意してもらわなくてもって感じですね。

割合ではなく要素数でデータがある場合は全体を合計で割って割合にすればよいだけですし。

以上で、ジニ不純度の説明を終えます。エントロピーと比べて一長一短あるのでお好みの方を使ってみてください。

シャノンの補助定理とその応用

今回も引き続き情報量(エントロピー)関係の記事です。直近、エントロピーに関する記事を複数書いていますが、その中でいくつか証明をしていない性質があります。シャノンの補助定理という面白い定理を使うとそれらが証明できるので見ていきます。

情報量の話なのでこの記事では底を省略して書いた対数関数の底は$2$とします。$\log=log_2$です。一方で、自然対数も使うのでそれは$\ln=\log_e$と書きます。

早速、シャノンの補助定理を見ていきます。

定理

二つの確率事象系$A = \{a_k | k=1,\dots, n \}$と$B = \{b_k | k=1,\dots, n\}$の確率をそれぞれ$p_k$, $q_k$とします。つまり、
$\sum_{k=1}^n p_k = 1$, $\sum_{k=1}^n q_k = 1$です。この時次の関係が成り立ちます。

$$ – \sum_{k=1}^n p_k\log{p_k} \leq – \sum_{k=1}^n p_k\log{q_k}.$$

以上が、シャノンの補助定理の主張です。左辺は$A$のエントロピーで、右辺は対数関数の中身だけ別の確率事象系の確率になっていますね。

証明

定理を証明する前に、自然体数に関する$\ln{x} \leq x – 1$ $(x > 0)$という関係を使うのでこれを先に証明します。

正の実数$x$に対して、$f(x) = x -1 – \ln{x}$ とおきます。すると$f'(x) = 1 – \frac{1}{x}$となります。$f'(x)$の値を見ていくと、 $0 < x < 1$ の時$f'(x)<0$であり、$f'(1)=0$、そして$1 < x$の時$f'(x)>0$となるので、$x=1$で最小値をとり、その値は$f(1) = 1-1-\ln{1} = 0$です。

よって、$f(x) \leq 0$であり、$\ln{x} \leq x – 1$ が示されました。

ここから定理の証明に入ります。定理の左辺から右辺を引いたものを考えて、底の変換をし、先ほどの関係式を適用します。

$$ \begin{align}
– \sum_k p_k \log{p_k} + \sum_k p_k\log{q_k} &= \sum_k p_k \log{\frac{q_k}{p_k}}\\
&=\frac{1}{\ln{2}} \sum_k p_k \ln{\frac{q_k}{p_k}}\\
&\leq \frac{1}{\ln{2}} \sum_k p_k \left(\frac{q_k}{p_k} – 1 \right)\\
&= \frac{1}{\ln{2}} \sum_k (q_k – p_k)\\
&= \frac{1}{\ln{2}} \sum_k q_k – \frac{1}{\ln{2}} \sum_k p_k\\
&=0\end{align}.$$

よって、$ – \sum_{k=1}^n p_k\log{p_k} \leq – \sum_{k=1}^n p_k\log{q_k}$ が示されました。等号が成立するのは$p_k = q_k$の場合です。

この補助定理を使って、エントロピー関連の性質を見ていきます。

応用1. エントロピーの最大値

最初に見ていくのは、「平均情報量(エントロピー)の定義について」で見た、エントロピーの範囲、$0 \le H(a) \le \log{n}$ の最大値の方です。これは、$q_k = \frac{1}{n}$としてシャノンの補助定理を使うと証明できます。

$$H(a) = – \sum_k p_k\log{p_k} \leq – \sum_k p_k\log{\frac{1}{n}} = \log{n}$$

となります。

応用2. 相互情報量が0以上であること

次に証明するのは、「相互情報量について」で最後に書いた、相互情報量が非負の値を取るという話です。

これは、シャノンの補助定理をそのまま2次元の確率分布に拡張して使います。
$k = 1, \dots, n$, $l = 1, \dots, m$とし、$\sum_{k, l} p_{kl} = 1$、$\sum_{k, l} q_{kl} = 1$とすると、

$$- \sum_k \sum_l p_{kl} \log{p_{kl}} \leq – \sum_k \sum_l p_{kl} \log{q_{kl}}$$

ですから、

$$\sum_k \sum_l p_{kl} \log{\frac{p_{kl}}{q_{kl}}} \geq 0$$

となります。ここで、$p_{kl} = P(a_k, b_l)$とし、$q_{kl} = P(a_k) P(b_l)$とすると、

$$\sum_k \sum_l P(a_k, b_l) \log{\frac{P(a_k, b_l)}{P(a_k) P(b_l)}} \geq 0$$

となります。この左辺が相互情報量$I(A; B)$の定義そのものなのでこれが非負であることが示されました。

応用3. 条件付きエントロピーがエントロピーより小さいこと

3つ目はさっきの相互情報量の話から続けて導けるものです。「結合エントロピーと条件付きエントロピー」で、重要な性質として$H(A|B) \le H(A)$が成り立つと証明無しで言いました。

これはシンプルに、$I(A; B) = H(A) – H(A|B)$であり、$I(A; B) \geq 0$なので、移項するだけで、$$H(A) \geq H(A|B)$$ が言えます。

まとめ

この記事ではシャノンの補助定理とその応用をいくつか紹介しました。応用1,2,3で証明したことはそれぞれの記事で証明無しで紹介していたのを個人的にモヤモヤしていたので、それを解消できてよかったです。

エントロピーの最大値についても、相互情報量の非負性についても結構シンプルな主張なのですが、シャノンの補助定理無しで証明するのは結構難しいのでこの定理の偉大さを感じます。

ここしばらくエントロピーの理論面の記事が続いていますが、次あたりでPythonで実際に計算する方法を紹介してこのシリーズを完結させようかなと思っています。

相互情報量について

相互情報量の定義と意味

エントロピー(情報量)関係の記事の3つ目です。前回の記事の最後で、エントロピー、結合エントロピー、条件付きエントロピーの間に成り立つ関係式で4つの量が等しくなるという話をしました。
参考: 結合エントロピーと条件付きエントロピー

その値のことを、相互情報量と呼び、$I(A; B)$と書きます。($A$と$B$は対等で交換も可能なのに$I(A, B)$ではなく$I(A; B)$と書く理由が気になりますね。)

改めて式を書くとこうなります。

$$\begin{align} I(A;B) &= H(A)-H(A|B)\\
&=H(B)-H(B|A)\\
&=H(A)+H(B)-H(A, B)\\
&=H(A, B)-H(A|B)-H(B|A)\end{align}$$

現実的には、一番上の行の等式で定義することが多いようです。

事象$A$の不確実性である、$H(A)$から、事前情報として$B$を知っている場合の$A$の不確実性$H(A|B)$を引いているわけですから、$I(A;B)$は、事前情報として$B$を知ったことによって減った$A$の不確実性と考えられます。$B$について知ることで分かった$A$に関する情報の方がわかりやすいかもですね。

個人的には$H(A) = I(A; B) + H(A|B)$ と移項して、{Aの情報量} = {Bを知った時点で得られるAの情報量} + {Bを知った後にAを観察して初めてわかるAの情報量} と考える方が理解しやすいと思っています。

この相互情報量は情報理論において非常に重要な概念です。これは、2つの確率変数間の情報の共有度合いを測る指標として使用されています。

相互情報量の確率表現

実用上の話として、$H(A)-H(A|B)$のままでは計算しにくいので確率$P$の式として相互情報量を表すことを目指します。一番上の定義をスタートとして計算してみましょう。

$$\begin{align}I(A;B) &= H(A)-H(A|B)\\
&= -\sum_{A}P(a)\log{P(a)} – \left\{ – \sum_{A}\sum_{B} P(a, b) \log{P(a|b)} \right\}\\
&= -\sum_{A}\sum_{B}P(a, b)\log{P(a)} + \sum_{A}\sum_{B} P(a, b) \log{\frac{P(a,b)}{P(b)}} \\
&= \sum_{A}\sum_{B} P(a, b) \log{\frac{P(a, b)}{P(a)P(b)}}.
\end{align}$$

なかなか綺麗な形に導けましたね。この形で見ると$A$, $B$を入れ替えても同じ値になることなどが容易にわかります。

Wikipedia などもそうですが、この確率表現の形の方を相互情報量の定義として、先に挙げた等式類を定理として導く流儀もあるようです。というより僕自身も元々そちらで理解していました。

相互情報量の性質

相互情報量は常に非負の値を取ります。

$$I(A;B)\ge 0.$$

また、$A$と$B$が独立の時、$P(a, b) = P(a)P(b)$ ですから、$\log$の中身が$1$になり、値が$0$になるため、$I(A;B)=0$となります。

これに関してはここ最近の証明してない他の性質とまとめて別記事で証明を紹介します。

逆に、$I(A; B) \le H(A)$ と $I(A; B) \le H(B)$ もそれぞれ成立します。

$A$の例を見ていきますが、等号が成立するのは $H(A|B)=0$となる場合です。これは、$B$の情報を得た時点で$A$のことが完全にわかってAを知っても情報が増えないケースが該当します。

結合エントロピーと条件付きエントロピー

前回の記事に続いて平均情報量(エントロピー)の話です。

参考: 平均情報量(エントロピー)の定義について

今回の記事ではその2次元版とも言える2種類のエントロピーを紹介します。この2つについてはどうも情報量よりエントロピーと呼ぶことが多いようなのでそのように書きます。また、数式中に和の記号がたくさん登場し、添え字書いていくと大変だし、かえって混乱もまねくので、$\sum_{k=1}^n f(a_k)$ を $\sum_{A}f(a)$のようにシンプルに書きます。

そのため、$A = \{a_k\}$のエントロピーは次のように書きます。

$$H(A) = \sum_{A}P(a)\log{P(a)}.$$

参考にした本は前回の記事同様、平田先生の「情報理論のエッセンス」です。

結合エントロピー

まず紹介するのはエントロピーの2次元版ともいえる、結合エントロピーです。これは二つの事象の集合 $A={a_k}$ と $B={b_l}$に対して定義されます。$a_k$と$b_l$が発生する確率を$P(a_k, b_l)$ とすると、結合エントロピーは次のように定義されます。

$$H(A, B) = – \sum_{A}\sum_{B}P(a, b)\log{P(a, b)}.$$

ここで、$\sum_{A}\sum_{B}P(a, b)=1$であることに注意してください。

これはもう、$A$と$B$の組みあわせとして発生しうる事象を全て列挙してそれに対してエントロピーを計算したようなもので、とてもわかりやすい定義だと思います。

$A$と$B$の結果を知った時に得られると期待できる情報量、とも言えます。

条件付きエントロピー

次に紹介するのは条件付きエントロピーです。これは先ほどの結合エントロピーと少し状況が異なり、あらかじめ$B$に関する情報を知った上で後から$A$についての結果を得たらどれだけの情報量を得られるかを表すものです。

こちらは結合エントロピーと違って、定義をバシッと書いてもわかりにくいので例を取り上げて順番に説明します。

本の例では、$A$がサイコロの目で$\{1, 2, 3, 4, 5, 6\}$、$B$が偶数/奇数というものを取り上げていましたのでここでもそれにならいます。$b_1$が偶数で$b_2$が奇数です。

まず$A$のエントロピーは$H(A)=-\sum_{k=1}^{6}\frac{1}{6}\log{\frac{1}{6}}=\log{6}\fallingdotseq 2.58$ です。

ここで、サイコロをの目を見る前になぜか偶数か奇数かだけわかり、$b_1$ = 偶数だったとしましょう。その後に目を確認するとすると、$2, 4, 6$が出ている確率は$1/3$で、$1, 3, 5$が出ている確率は$0$ということになります。

つまりサイコロの目を見ることで得られる情報量は、

$H(A|b_1) = -\sum_{A}P(a|b_1)\log{P(a|b_1)} = -\sum_{a\in\{2, 4,6\}}\frac{1}{3}\log{\frac{1}{3}} = \log{3} \fallingdotseq 1.58$

となります。情報量が1だけ減っていますね。この1は、$\log{2}$相当の値で、事前に偶数か奇数かという情報を得たことで起こりうる確実性が$1/2$になったことを示しています。

同様に計算すると、事前に奇数だっと知らされた後に出た目を確認することによって得られる情報量$H(A|b_2)$も$\log{3}$だとわかります。

そして、いよいよ本題なのですが、ここで算出した$H(A|b_1)$, $H(A|b_2)$に、$B$の結果が$b_1$であるか$b_2$であるか、要するに偶数であるか奇数であるかの確率をかけて期待値を取ると、偶数か奇数かの事前情報がある場合の$A$のエントロピーの期待値$H(A|B)$を次のように求められます。

$$\begin{align}
H(A|B) &= H(A|b_1)P(b_1) + H(A|b_2)P(b_2)\\
&= \frac{1}{2}\log{3} + \frac{1}{2}\log{3} = \log{3} \fallingdotseq 1.58
\end{align}$$

この$H(A|B)$が条件付きエントロピーです。もっと一般の事例について定義し、$P$の式として表記できるよう変形しておくと、次のようになります。

$$\begin{align}
H(A|B) &= \sum_{B} H(A|b)P(b)\\
&= -\sum_{A}\sum_{B} P(a|b)P(b)\log{P(a|b)}\\
&= -\sum_{A}\sum_{B} P(a, b)\log{P(a|b)}.\\
\end{align}$$

結合エントロピーとぱっと見似ていますがよく見ると違いますね。$\log$の中身が結合確率ではなく条件付き確率になっています。

また、重要な性質として、次の大小関係が成り立ちます。

$$H(A|B) \le H(A).$$

各種エントロピー間の関係式

前回の記事も含めて3種類のエントロピーを定義したので、それぞれの間にある関係式を見ていきます。まず紹介するのは $H(A, B) = H(B) + H(A|B)$です。情報理論のエッセンスだと、ベン図みたいなのを描いてイメージで説明されていますが、計算すると次ようになります。

まず定義から次のようになります。

$$H(A, B) = – \sum_{A}\sum_{B}P(a, b)\log{P(a, b)}.$$

ここで、$P(a, b) = P(a | b) P(b)$であることと、$\log$内の積は和に分けれるので、次のようになります。

$$H(A, B) = – \sum_{A}\sum_{B}P(a, b)\log{P(a|b)} – \sum_{A}\sum_{B}P(a, b)\log{P(b)}$$

前半は$H(A|B)$の定義が出てきましたね。後半は、$\sum_{A}P(a, b) = P(b)$であることに注意すると、次のようになります。

$$H(A, B)= H(A| B) – \sum_{B}P(b)\log{P(b)} = H(A|B) + H(B)$$

$a, b$ 入れ替えて、$P(a, b) = P(b | a) P(a)$を使うと、全く同じようにして、$H(A, B) = H(A) + H(B|A)$ も示せます。

以上で、

$$H(A, B) = H(A) + H(B|A) = H(B) + H(A|B)$$

が示せました。これの中辺と右辺に着目して移行すると次も言えます。

$$H(A) – H(A|B) = H(B) – H(B|A).$$

また、左辺と右辺に着目し、$H(A, B)$と$H(A|B)$をそれぞれ移行して両辺に$H(A)$を足すとつ次の式が出てきます。

$$H(A) – H(A|B) = H(A) + H(B) – H(A, B).$$

左辺と中辺に着目して、$H(B|A)$を移行して両辺から$H(A|B)$を引くと、次の式が出てきます。
$$H(A) – H(A|B) = H(A, B) – H(A|B) – H(B|A).$$

以上をまとめると、次の4つの値は違いに等しいと言えます。
$H(A) – H(A|B)$
$H(B) – H(B|A)$
$H(A) + H(B) – H(A, B)$
$H(A, B) – H(A|B) – H(B|A)$

これを相互情報量と呼ぶのですが、長くなってきたので次の記事で詳しく説明します。

まとめ

今回の記事では結合エントロピーと条件付きエントロピーを紹介しました。

結合エントロピーの方は多次元へのシンプルな拡張というイメージでしたが、条件付きエントロピーはすでに何かしらの情報を得てる状況で新たに何かを知ったらどれほどの情報を得られるかを示すものでありなかなか重要な概念です。

また、結合エントロピーの方も重要ではないというわけではありません。実際にプログラム等で値を計算をするときは、条件付きエントロピーを直接算出する関数が主要ライブラリで提供されいていないという事情により、頻繁に算出することになります。

平均情報量(エントロピー)の定義について

はじめに

データサイエンスや機械学習関連の本を読んでいると時々エントロピーという概念に出逢います。例えば、決定木の分類の品質の測定などはそうですね。scikit-learnでは引数でジニ不純度とエントロピー、ログロスの中から選べるので必須でも無いのですが。

個人的な話ですが、このエントロピーは定義がいまいちしっくりこなかったのであまり使ってきませんでした。scikit-learnの決定木もcriterionのデフォルトは”gini”だし、こちらの方が扱いやすく感じていたというのもあります。

ただ、最近平田先生の「情報理論のエッセンス」という本を読みましてこのエントロピーについてだいぶ理解が深まったので軽くまとめておきます。

平均情報量の定義

この記事では離散な事象の場合を取り扱います。最初に平均情報量(エントロピー)の定義を書いておきます。

事象の集合$A=\{a_k\,|\,k=1, 2, \dots, n \}$の各事象$a_k$に対して、その事象が発生する確率を$P(a_k)$とすると、平均情報量$H(A)$は次のように定義されます。

$$H(A) = -\sum_{k=1}^nP(a_k)\log{P(a_k)}.$$

この平均情報量は、式の形が熱力学等のエントロピーと同じなので、情報理論においてもエントロピーと呼ばれています。

この記事では、平均情報量がなぜこのように定義されるのかというのを見ていきます。

熱力学の便利な式を適当に持ってきたわけでもなく、何となく$\log$を使っているわけでもなくきちんと背景があってこのように定義されたのだ、というのが伝われば幸いです。

情報量の数値化について

そもそものモチベーションは情報を数値化することにあります。何かの事象を観測したり、ニュースを聞いたりした時にその結果や内容がどれだけ「予想外」もしくは「驚き」であったかを数値化したものです。

先述の情報理論のエッセンスでは、次の二つのニュースが例として挙げられていました。
(a) 東京に雪が降りました。
(b) 北極に雪が降りました。

(b)の方はありふれた話であるのに対して、(a)の方は珍しいことであって(a)のニュースの方がニュース価値が大きいものであると考えられます。

このニュースとしての価値を定量的に定義しようというのが、そもそもの情報量の目的です。

情報量が満たすべき性質とそれを満たす式

ここから、具体的に事象$x$の情報量$i(x)$について考えていきますが、そのためにこれがどんな性質を持つべきなのか見ていきます。

まず考えておかないといけないのは、受け取り手が得られる情報としての価値にはその人の主観が入り得ることです。東京の天気には一切興味がなく、北極にはすごく関心があるという人にとっては(b)のニュースの方が価値があるかも知れません。これは「主観的立場」による情報の評価ですが、情報理論で扱う情報量はそうではなく、誰でも同じように評価できるようにするために「客観的立場」で情報を評価します。

客観的に評価するために、情報量を定義する際はその内容ではなく、そのニュースが発生する確率を元に価値を評価します。そのため、「情報量は確率の関数である」というのが一つ目の条件です。つまり次のように書けます。

$$i(x) = i(P(x)).$$

そして、起こりにくい事象ほどそれについて知った時の情報価値は大きいということから、$P(a)<P(b)$ならば$i(a) > i(b)$となっていることが求められます。つまり、二つ目の条件は次のようにいえます。

$i(P(x))$は$P(x)$の減少関数である。

この時点では減少関数なんていくらでもあるので情報量がどうとでも定義できてしまいます。

そこで最後に、情報の加法性という性質を要求します。これは独立した二つの情報を知った時に、その情報を同時に知った場合と、順番に知った場合、またその逆順に知った場合で得られる情報の量は同じである、ということです。

例えば、トランプを1枚引いて、以下のように確認した場合に最終的に得られる情報は全部一緒ということです。
– そのカードはスペードの10だった。
– そのカードの数字は10だった。そしてマークはスペードだった。
– そのカードのマークはスペードだった。そして数字は10だった。

スペードだったという情報を$a$、 数字が10だったという情報を$b$とすると、
$i(a\land b) = i(a) + i(b) = i(b) + i(a)$ であって欲しいということです。二つの情報が独立である場合、$P(a, b) = P(a) P(b)$ですから、この性質は次のように書くことができ、これが三つ目の条件です。

$$i(P(x)P(y)) = i(P(x)) + i(P(y)).$$

これを満足する関数を考えると、対数関数しかなくなります。

先に出てる条件で減少関数としないといけませんから、底を$s$として次のようになります。

$$i(x) = i(P(x)) = -\log_s{P(x)}$$

底の$s$としては$2$,$e$,$10$など考えられますが、情報理論においては$s=2$が使われます。ここから底は省略して$2$とします。

コインの裏表やスイッチのON/OFFなど、$1/2$の確率で得られる情報の情報量がちょうど$1$になります。また、6面のサイコロを転がして何か目が出たら、そこから得られる情報量は$\log{6}$となります。

平均情報量の定義

ここまでで、何かニュースを得たり事象を観測したらどの程度の情報が得られるかを定量化することができました。

次に、そのニュースを得る前に、そのニュースを得たらどの程度の情報量が得られるか、の期待値をを考えます。つまり、先述の例で言えば、「東京に雪が降るかどうか」を確認したらどの程度の情報量が得られるのかの期待値を定量化しようという話です。

これは、各事象ごとにその事象が発生する確率と、得られる情報の積をとって足し合わせることで得られます。通常の期待値の定義ですね。つまり平均情報量は次のようになります。

$$\begin{align}H(A) &= \sum_{k=1}^n i(a_k)P(a_k)\\
&=-\sum_{k=1}^n P(a_k)\log {P(a_k)}
\end{align}$$

これでこの記事の最初にあげた定義式が出てきました。

元々$-p\log{p}$ってなんだ?と思ってたのですが、導出を追いかけてみるとこれしか無いって感じる自然な定義ですね。

平均情報量の性質

基本的な性質を一つ紹介しておきます。平均情報量(エントロピー)$H(A)$は次の関係を満たします。

$$0 \le H(a) \le \log{n}.$$

左側の等号は $^{\exists}l \; P(a_l) = 1, k\neq l \Rightarrow P(a_k)=0$ の時に成立します。1種類の結果しか生じ得ない実験をやっても得られる情報量は0であるってことですね。

右側の統合は、$^\forall k P(a_k)=\frac{1}{n}$の時に成立します。これは発生しうる$n$種類の結果が全て等しい確率で起こりうる場合ですね。

平均情報量のもう一つの見方

特にエントロピーと呼ばれる時はそうなのですが、平均情報量は「不確実性の尺度」として理解されることも多くあります。熱力学ではそうですね。

不確実性と得られる情報量では全く逆のものに見えるのですが、これが同一の式で定義されているのは、「得られる情報量」を「不確実性の減少量」と考えられるからです。

ある実験をする前に、この事象にはこれだけ不確実性があるぞ、と評価してるのが熱力学等のエントロピーで、事象を実際に観測してこれだけの情報が得られたと言っているのが情報理論のエントロピーって感じですかね。

まとめ

ここまでで、情報の量を定式化したいというモチベーションから、定式化するとしたらどんな性質を持たないといけないか、そしてどのような式で定義できるかを見てきました。

平均情報量(エントロピー)の定義にについての話は以上になります。次回以降の記事でもう少し発展させて条件付きエントロピーや相互情報量について取り扱っていこうと思います。

2標本コルモゴロフ–スミルノフ検定を試す

前回の記事の続きです。
参考: 1標本コルモゴロフ–スミルノフ検定について理解する

コルモゴロフ-スミルノフ検定(K-S検定)は、一つの標本が何かしらの確率分布に従っているかどうかだけではなく、二つの標本について、それらが同一の確率分布から生成されたものかどうかの判定にも利用可能です。

Wikipediaを見ると、KS検定は1標本より2標本の時の利用の方が推されてますね。2標本それぞれの確率分布が不明となるとノンパラメトリックな手法が有効になるのでしょう。

前の記事みたいに数式をつらつらと書いていこうかと思ったのですが、1標本との違いは、統計量を計算するときに、経験分布と検定対象の累積分布関数の差分の上限(sup)を取るのではなく、その二つの標本それぞれの経験分布に対して、差分の上限の上限(sup)を取るというだけです。KS検定について調べるような人でこの説明でわからないという人もいないと思うので、数式等は省略します。(前回の記事を参照してください。)

便利な特徴としては、2標本のそれぞれのサイズは等しくなくても使える点ですね。
ただ、色々試した結果、それなりにサンプルサイズが大きくないとあまり使い物にならなさそうです。

この記事では、Pythonで実際に検定をやってみます。
前回のt分布からの標本と正規分布は、標準偏差が違うとはいえ流石に似すぎだったので、今回はちょっと変えて、カイ2乗分布と、正規分布を使います。
カイ2乗分布の自由度は4とし、すると期待値4、分散8になるので、正規分布の方もそれに揃えます。期待値や分散が違うならt検定やF検定で十分ですからね。

とりあえず分布関数を作り、可視化して見比べておきます。

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


chi2_frozen = chi2.freeze(df=4)  # 自由度4のカイ2乗分布 (期待値4, 分散8)
norm_frozen = norm.freeze(loc=4, scale=8**0.5)  # 正規分布 (期待値4, 分散8)

# 確率密度関数を可視化
xticks = np.linspace(-10, 20, 301)
fig = plt.figure(facecolor="w")
ax = fig.add_subplot(1, 1, 1)
ax.plot(xticks, chi2_frozen.pdf(xticks), label="カイ2乗分布")
ax.plot(xticks, norm_frozen.pdf(xticks), label="正規分布")
ax.legend()
plt.show()

出力はこちら。これはちゃんと検定で見分けて欲しいですね。

続いて、それぞれの分布から標本を取ります。サンプルサイズが違っても使える手法なので、意図的にサンプルサイズ変えました。

# それぞれの分布から標本を生成
chi2_samples = chi2_frozen.rvs(size=200, random_state=0)  # カイ2乗分布からの標本
norm_samples = norm_frozen.rvs(size=150, random_state=0)  # 正規分布からの標本

さて、これでデータが取れたので、2標本KS検定やっていきます。
非常に簡単で、ks_2samp ってメソッドに渡すだけです。
ドキュメント: scipy.stats.ks_2samp — SciPy v1.9.0 Manual

例によって、alternative 引数で両側検定/片側検定などを指定できますが、今回両側でやるので何も指定せずにデフォルトのまま使います。

帰無仮説は「二つの標本が従う確率分布は等しい」です。有意水準は0.05とします。

結果がこちら。

print(ks_2samp(chi2_samples, norm_samples))
# KstestResult(statistic=0.17, pvalue=0.012637307799274388)

統計量とp値が表示されました。p値が0.05を下回りましたので、帰無仮説を棄却し、二つの標本が従う確率分布は異なっていると判断します。

サンプルサイズさえそこそこ確保できれば、非常に手軽に使えて便利な手法だと思いました。

あとは個人的には、KS検定って5段階や10段階評価のような離散分布でも使えるのかどうか気になっています。まぁ、5段階の場合は単純にカイ2乗検定しても良いのですが、10段階以上になってくると自由度が高くてカイ2乗検定があまり使いやすくないので。何か情報がないか追加で探したり、自分でシミュレーションしたりしてこの辺をもう少し調べようと思います。

1標本コルモゴロフ–スミルノフ検定について理解する

何らかのデータ(標本)が、特定の確率分布に従ってるかどうかを検定したいことって頻繁にあると思います。そのような場合に使えるコルモゴロフ–スミルノフ検定(Kolmogorov–Smirnov test, KS検定)という手法があるのでそれを紹介します。

取り上げられている書籍を探したのですが、手元に見当たらなかったので説明はWikipediaを参照しました。
参考: コルモゴロフ–スミルノフ検定 – Wikipedia

1標本と、特定の確率分布についてその標本が対象の確率分布に従っていることを帰無仮説として検定を行います。

この記事では、scipyを使った実装と、それを理解するための検証をやっていきたいと思います。

とりあえずこの記事で使うライブラリをインポートし、データを準備します。
標本は、自由度5のt分布からサンプリングし、それが正規分布に従ってないことを検定で見ていきましょう。scipyのパラメータはどちらの分布もloc=10, scale=10 としました。

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


norm_frozen = norm.freeze(loc=10, scale=10)  # 検定する分布
t_frozen = t.freeze(df=5, loc=10, scale=10)  # 標本をサンプリングする分布
t_samples = t_frozen.rvs(size=1000, random_state=0)  # 標本を作成

# 各データを可視化
xticks = np.linspace(-40, 60, 1001)
fig = plt.figure(facecolor="w")
ax = fig.add_subplot(1, 1, 1)
ax.plot(xticks, norm_frozen.pdf(xticks), label="正規分布")
ax.plot(xticks, t_frozen.pdf(xticks), label="t分布")
ax.hist(t_samples, density=True, bins=50, alpha=0.5, label="t分布からの標本")
ax.legend()
plt.show()

作成された図がこちらです。流石にt分布と正規分布は似ていますね。ぱっと見でこの標本が正規分布に従ってないと言い切るのは難しそうです。

それでは、早速scipyで(1標本)KS検定をやっていきましょう。これは専用の関数が用意されているので非常に簡単です。
参考: scipy.stats.kstest — SciPy v1.9.0 Manual

rvs引数に標本データ、cdf引数に検定したい分布の累積分布関数を指定してあげれば大丈夫です。alternative で両側検定、片側検定などを指定できます。今回は両側で行きます。有意水準は0.05としましょう。

ks_result = kstest(rvs=t_sample, cdf=norm_frozen.cdf, alternative="two-sided")

print(ks_result)
# KstestResult(statistic=0.04361195130340301, pvalue=0.04325168699194537)

# 統計量と引数を変数に入れておく
ks_value = ks_result.statistic
p_value = ks_result.pvalue

はい、p値が約0.043 で0.05を下回ったので、この標本が正規分布にし違うという帰無仮説は棄却されましたね。

ちなみに、標本をサンプリングした元のt分布でやると棄却されません。これも想定通りの挙動ですね。

print( kstest(rvs=t_sample, cdf=t_frozen.cdf, alternative="two-sided"))
# KstestResult(statistic=0.019824350810698554, pvalue=0.8191867386190135)

一点注意ですが、どのような仮説検定でもそうである通り、帰無仮説が棄却されなかったからと言って正しいことが証明されたわけではありません。(帰無仮説は採択されない。)
特にKS検定については、標本サイズを変えて何度も試してみたのですが、標本サイズが小さい時は誤った帰無仮説を棄却できないことがかなりあるようです。

さて、scipyで実行してみて、この検定が使えるようになったのですが、より理解を深めるために、自分で統計量を計算してみようと思います。

KS検定の統計量の計算には、経験分布というものを使います。要するに標本から作成した累積分布関数ですね。数式として書くと標本$y_1, y_2,\dots, y_n$に対する経験分布$F_n$は次のようになります。

$$F_n(x) = \frac{\#\{1\le i\le n | y_i \le n\}}{n}$$

要するに$x$以下の標本を数えて標本のサイズ$n$で割ってるだけです。

そして、検定したい分布の分布関数$F$とこの$F_n$に対して、KS検定の統計量は次のように計算されます。

$$\begin{align}
D^{+}_n &= \sup_x\left(F_n(x)-F(x)\right)\\
D^{-}_n &= \sup_x\left(F(x)-F_n(x)\right)\\
D_n &= \max(D^{+}_n, D^{-}_n)
\end{align}$$

max(最大値)ではなく、sup(上限)で定義されているのがポイントですね。
とはいえ、分布関数の方が一般的な十分滑らかな関数であれば、$D^{+}_n$の方はmaxでもsupでも同じですが,$D^{-}_n$の方はsupが効いてきます。

この$D_n$の直感的なイメージを説明すると、$F$と$F_n$の差が一番大きくなるところの値を取ってくることになるでしょうか。

経験分布の方を実装して可視化してみます。

# 経験分布関数
def empirical_distribution(x, samples):
    return len([y for y in samples if y <= x])/len(samples)


xticks_2 = np.linspace(-40, 60, 100001)  # メモリを細かめに取る
# 経験分布関数の値
empirical_cdf = np.array([empirical_distribution(x, t_samples) for x in xticks_2])

# 可視化
fig = plt.figure(facecolor="w")
ax = fig.add_subplot(1, 1, 1)
ax.plot(xticks_2, norm_frozen.cdf(xticks_2), label="正規分布")
ax.plot(xticks_2, empirical_cdf, label="経験分布")
ax.legend()
plt.show()

出てきた図がこちらです。標本サイズが大きいのでよくみないとわかりませんが、オレンジの方は階段状の関数になっています。

この、青線(検定する正規分布)とオレンジの線(標本からの経験分布)が一番離れたところの差を統計量にしましょう、というのが考え方です。

では$D^+_n$から計算しましょう。実は(今回の例では)青線が滑らかなのに対して、オレンジの線が階段状になっているので、この段が上がるポイント、つまり標本が存在した点での値だけ調べると$D^+_n$が計算できます。

dplus = max([empirical_distribution(x, t_samples) - norm_frozen.cdf(x) for x in t_samples])
print(dplus)
# 0.031236203108694176

次に$D^-_n$の方ですが、これは少し厄介です。サンプルが存在する点ではなく、その近傍を調べて階段の根元の値と累積分布関数の差を取りたいんですよね(ここでmaxではくsupが採用されていた意味が出てくる)。1e-10くらいの極小の定数を使って計算することも考えましたが、経験分布関数の方を少しいじって計算してみることにしました。どういうことかというと、x以下の要素の割合ではなくx未満の要素の割合を返すようにします。これを使うとサンプルが存在した点については階段の根元の値が取れるようになるので、これを使って$D^+_n$と同様に計算してみます。

def empirical_distribution_2(x, samples):
    return len([y for y in samples if y < x])/len(samples)


dminus = max([norm_frozen.cdf(x) - empirical_distribution_2(x, t_samples) for x in t_samples])
print(dminus)
# 0.04361195130340301

さて、必要だった統計量は$D^+_n$,$D^-_n$の最大値です。これをライブラリが計算してくれた統計量と比較してみましょう。一致しますね。

print("自分で計算したKS値:", max([dplus, dminus]))
# 自分で計算したKS値: 0.04361195130340301
print("scipyのKS値:", ks_vakue)
# scipyのKS値: 0.04361195130340301

さて、統計量がわかったら次はp値を計算します。ウィキペディアを見ると、この統計量は無限和を含むそこそこ厄介な確率分布に従うことが知られているそうです。
これからp値を求めるのは流石に手間なのと、幸い、scipyに実装があるので(scipyの検証にscipyを使うのは不本意ですが)ここは甘えようと思います。

非常にありがたいことに、kstwoksoneという両側検定、片側検定に対応してそれぞれ分布関数を用意してくれています。

統計量は求まっていますし、分布関数もあるのでp値はすぐ出せます。

print(kstwo.sf(max([dplus, dminus]), n=1000))
# 0.04325168699194537

# 参考 kstestから取得したp値
print(p_value)
# 0.04325168699194537

最後、少し妥協しましたが、追々kstwo分布についても自分でスクラッチ実装して検証しておこうと思います。

今回の検証でコルモゴロフ–スミルノフ検定についてだいぶ理解が深まったので、これからバシバシ使っていこうと思います。

1標本だけでなく2標本でそれぞれの分布が等しいかどうかを検定するって使い方もできるので、次の記事はそれを取り上げようかなと思っています。