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

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

共役事前分布とは

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

例えば、ベータ分布は二項分布の共役事前分布なので、二項分布のパラメーター$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$のカイ二乗分布になります。
カイ二乗分布もこのブログで紹介しているので確率密度関数の式の形を見比べてみてください。
参考: カイ二乗分布について

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

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

最高密度区間 HDI (Highest Density Interval) について

これはpyMCに限った話ではないのですが、pyMC + ArviZ を使っていると頻繁に目にする指標なのでほぼpyMC関連の記事6記事目と見ていただいて良いと思います。

pyMCの結果を見ているとhdiという指標を頻繁に目にします。僕はこれについてあまり馴染みがなかったので説明を残しておこうと思います。

各記事のpyMCのサンプリング結果のsummaryをDataFrameにしたときも登場していますし、可視化したらグラフ中にも登場しますね。

参考: ArviZを使ってpyMCの推論結果を可視化する

このHDIは、確率分布の区間を切り取ったものです。指定の確率、(一般的には95%だそうですが、) ArviZでは94%がデフォルトで指定されており、その区間に含まれる確率が94%となるような区間を示します。

似たような概念で信頼区間や予測区間というのもありますが、違いはその区間の切り取り方です。信頼区間や予測区間では大抵、区間を切り取るとき、確率分布の上側と下側から同じ確率を切り抜きます。94%区間を取りたいのであれば、下から3%と上から3%を外して残りを取り出します。

それに対して、最高密度区間(HDI)では、その名の通り、確率密度が高い部分から優先的に確保して、指定の確率、今の例であれば94%となるように区間を切り出します。

山が一つあるような形の確率分布$f(x)$であれば、次を満たす区間$[a, b]$だと言えるでしょうか。$f(x)$が連続とか滑らかとか適切な過程が他にも必要かもしれません。

$$\begin{align}
&\int_{a}^{b} f(x)dx = 0.94\\
&f(a) = f(b)\\
&a < c < b \Rightarrow f(a) = f(b) < f(c)
\end{align}$$

注意として、確率分布が複数の山を持つような形だった場合、HDIは単一の区間ではなく、ふくすの区間に分かれて定義されることがあります。pyMCやArviZがそのような例に対応しているかどうかはちょっと確認しきれていません。(まだそういう例に出くわしたことがないので)

HDIの良い性質として、その指定の確率を含む区間の切り取り方としては一番狭い幅の区間になるというものがあります。ベイズ推論に限った話ではなく、何か数値を予想するのであればできるだけ狭いレンジに絞り込んで推論したいのでこれはありがたいですね。