ポアソン混合モデル　その1：ポアソン分布とガンマ分布、それらの共役性

※このページでは混合分布の理論は一切登場しません。※

ポアソン混合モデルでは、ポアソン混合分布、つまり、いくつかのポアソン分布を合わせた分布を扱う。

下記に、3つのポアソン分布を混ぜたポアソン混合分布を描いてみた。

import numpy as np
from scipy.special import factorial
import matplotlib.pyplot as plt

x = np.linspace(start=0, stop=50, num=1000)
# ポアソン分布の尤度関数だけ作っておいた
uni_poisson = lambda x, lambda_: (lambda_ ** x) / factorial(x) * np.exp(-lambda_)
poisson1 = uni_poisson(x, 3)
poisson2 = uni_poisson(x, 10)
poisson3 = uni_poisson(x, 40)

# 合計が1になるような3つの重み(2次元標準シンプレックス)を作って足し合わせる
mixed_poisson = poisson1 * 0.2 + poisson2 * 0.3 + poisson3 * 0.5

plt.plot(x,mixed_poisson)

下記ではポアソン混合モデルの基礎となるポアソン分布や関連するガンマ分布について実装を行いつつ、その性質を断片的に確認していく。

###ポアソン分布

ポアソン分布(Poisson distribution)は非負の整数を発生させる離散確率分布である。
尤度関数は

Poi(x|λ) = \frac{λ^x}{x!} e^{-λ}

である。ここで、$λ \in \mathbb{R}^+$ はポアソン分布のパラメータ。

ちなみに、$Poi(x|λ)$ の期待値は $λ$ である(証明略)。

さて、ベイズ推論で使いがちな対数尤度関数は

\ln Poi(x|λ) = x \ln λ - \ln x! - λ

である。

パラメータ $λ = 5, 10, 20$ についてグラフを書いてみる。

import numpy as np
from scipy.special import factorial

def Poisson(x, lambda_):
  return (lambda_ ** x) / factorial(x) * np.exp(-lambda_)

# 描画用関数
import matplotlib.pyplot as plt
def show_Poisson(label, lambda_):
  lambda_ = lambda_
  x = np.linspace(start=0, stop=30, num=1000)
  y = Poisson(x, lambda_)
  plt.title(label=label, loc='center')
  plt.plot(x,y)

# 描画するよ
plt.figure(figsize=(12,3))

plt.subplot(131)
show_Poisson(lambda_ = 5, label='lambda_ = 5')

plt.subplot(132)
show_Poisson(lambda_ = 10, label='lambda_ = 10')

plt.subplot(133)
show_Poisson(lambda_ = 20, label='lambda_ = 20')

$Poi(x|λ)$ の期待値が $λ$ ってのは確かに、という感じ。

ガンマ分布

ポアソン混合モデルの文脈でガンマ分布が登場してくるときは、「ポアソン分布のパラメータ $ λ $ はガンマ分布に従うと仮定する」ってのが定石。解析的に解く場合はね。

ガンマ分布(Gamma distribution)は正の実数を生成させる連続確率分布。尤度関数は

Gam(λ|a, b) = C_G(a, b) λ^{a - 1} e^{- b λ}

である。ここで、$a, b \in \mathbb{R}^+$ はガンマ分布のパラメータ。

ちなみに、$Gam(x|λ)$ の期待値は $ \frac{a}{b} $ である(証明略)。

また、 $ C_G(a, b) $はガンマ分布の正規化項であり、以下の式で表される。

C_G(a, b) = \frac{b^a}{Γ(a)}

ちなみに、ガンマ関数の定義と重要な式を紹介しておく。

Γ(x) = \int t^{x - 1} e^{-t} dt \\
Γ(x + 1) = x Γ(x) \\
Γ(1) = 1

ベイズ推論で使いがちな対数尤度関数は

$$ \ln Gam(λ|a, b) = (a - 1) \ln λ - b λ + \ln C_G(a, b)$$

である。

パラメータ $ (a,b) = (1, 0.5), (2, 0.5), (3, 0.5), (5, 1)$ についてグラフを書いてみる。

import numpy as np
from scipy.special import gamma

def Gamma(lambda_, a, b):
  main = (lambda_ ** (a - 1)) * (np.exp(-b * lambda_))
  c_g = (b ** a) / gamma(a)
  return main * c_g

import matplotlib.pyplot as plt

a_b = [(1, 0.5), (2, 0.5), (3, 0.5), (5, 1)]

# 描画用関数
def show_Gamma(a_b):
  a, b = a_b[0], a_b[1]
  x = np.linspace(start=0, stop=20, num=2000)
  y = Gamma(x, a, b)
  # (a,b)の組と、期待値('Exp')を表示
  plt.title(label=f'(a,b) = {a_b}, Exp = {a/b}', loc='center')
  plt.ylim([0, 0.5])
  plt.plot(x,y)

# 描画するよ
plt.figure(figsize=(12,3))

plt.subplot(141)
show_Gamma(a_b[0])

plt.subplot(142)
show_Gamma(a_b[1])

plt.subplot(143)
show_Gamma(a_b[2])

plt.subplot(144)
show_Gamma(a_b[3])

期待値が$ \frac{a}{b} $ っていうのは図からもなんとなくわかる。

ちなみに、ガンマ分布のパラメータ $ b $ の逆数 として $ θ = \frac{1}{b} $ が使われることもある。

詳細は一旦おいておくが、$ θ $ が使われる場合のガンマ分布の尤度関数は下記になる。

$$ Gam(λ|a, θ) = \frac{1}{θ^a Γ(a)} λ^{a - 1} e^{- \frac{λ}{θ}} $$

ポアソン分布とガンマ分布の共役性

まず、「共役性」についての説明。ここではポアソン分布で説明。

パラメータは $ λ $。

所与データを$ \mathcal{D} $ (例えば、一次元データ $ (x_1, x_2, ... , x_n) $) と 表す。

事前分布 $ p(λ) $ と事後分布 $ p(λ|\mathcal{D}) $ が同じ関数形になるとき、事前分布 $ p(λ) $ は、尤度関数 $ p(\mathcal{D}|ω) $ に対する共役事前分布と呼ぶ。

ポアソン分布を例にとって説明する。ここでは混合モデルでなく単一のポアソン分布を例にとる。

尤度関数：

 \begin{align}
p(\mathcal{D}|λ) 
&= p(x_1, x_2, ... , x_n|λ) \\
&= \prod_{n}^{N} Poi(x_n|λ) \\
&= \prod_{n}^{N} \frac{λ^{x_n}}{x_n!} e^{-λ} \\
\end{align}

パラメータ $ λ $ の事前分布：

\begin{align}
p(λ) 
&= Gam(λ|a, b) \\
&= C_G(a, b) λ^{a - 1} e^{- b λ} \\
\end{align}

ここまでは復習+α。

尤度関数の式では、すべての $ x_n $ はが同一のポアソン分布に従って生成され、パラメータ $ λ $が与えられたもとでは $ x_i $ と $ x_j $ は独立と仮定している。

この仮定により、
$$ p(\mathcal{D}|λ) = \prod_{n}^{N} Poi(x_n|λ) $$

と、尤度関数の積に分解できる。

このとき、事後分布は

\begin{align}
p(λ|\mathcal{D})
&= \frac{p(\mathcal{D}, λ)}{p(\mathcal{D})} \\
&= \frac{p(\mathcal{D}|λ) p(λ)}{p(\mathcal{D})} \\
&\propto p(\mathcal{D}|λ) p(λ) \\
&= p(\mathcal{D}|λ) ・ Gam(λ|a, b) \\
&= \prod_{n}^{N} Poi(x_n|λ) ・ Gam(λ|a, b) \\
&= \prod_{n}^{N} \frac{λ^{x_n}}{x_n!} e^{-λ} ・C_G(a, b) λ^{a - 1} e^{- b λ} \\
&\propto \prod_{n}^{N} λ^{x_n} e^{-λ} ・ λ^{a - 1} e^{- b λ} \\
&= λ^{\sum_{n}^{N}x_n + a - 1} e^{- (N + b)λ} \\
&\propto Gam(λ| \sum_{n}^{N}x_n + a - 1, N + b) \\
\end{align}

という感じで、事後分布がガンマ分布になった。

ここで、注意事項。

・事後分布はパラメータ $ λ $ の関数とみているため、$ λ $ に依らない定数は式から追い出している。

・パラメータ$ λ $ に関する項（ここでは $ λ^{なにか} $ と $ e^{なにか × λ} $）が決まった時点で、ガンマ分布の形は決まる。

・ガンマ分布の確率密度関数
$$ Gam(λ|a, b) = C_G(a, b) λ^{a - 1} e^{- b λ} $$
のうち、$ λ $ に依らない $ C_G(a, b) $ は分布の積分値を1にするための定数であり、パラメータ $ λ $ に関する項が決まった時点で $ C_G(a, b) $ の値も確定するからあまり考えなくても大丈夫。

・尤度関数の対数を表現を使うとシンプルで見やすいので、興味ある方はぜひ調べるか手計算で事後分布を求めてみてね。

ポアソン分布で予測する

データ $ \mathcal{D}(x_1, x_2, ... , x_n) $ の分布をポアソン分布でうまく表したいと考える。

上述の通り、各データ $ x_n $ について

尤度関数：
$$ p(x_n|λ) = Poi(x_n|λ) $$

パラメータ $ λ $ の事前分布：
$$ p(λ) = Gam(λ|a, b) $$

と仮定したとき、

パラメータ $ λ $ の事後分布：

$$ p(λ|\mathcal{D}) = Gam(λ| \sum_{n}^{N}x_n + a - 1, N + b) $$

である。

最尤推定みたいに一つの実現値に固定することなくパラメータの事後分布を求めることができた、というのはベイズビギナーの私としては興奮ポイント。

未知データ $ x_n $ についても分布を求めることができ、その方法を以下に示す。

\begin{align}
p(x_*) 
&= \int {p(x_*|λ) p(λ)} dλ\\
&= \int {Poi(x_*)・Gam(λ| \sum_{n}^{N}x_n + a - 1, N + b)} dλ
\end{align}

データに基づいてそれらしく更新されたガンマ分布を採用。見づらいので、

a' = \sum_{n}^{N}x_n + a - 1 \\
b' = N + b

と置き換えておこう。

\begin{align}
p(x_*) &= \int{p(x_*|λ) p(λ)} \\
&= Poi(x_*)・Gam(λ| a', b') \\
&= \int{\frac{λ^{x_*}}{x_*!}e^{-λ}・C_G(a',b')λ^{a'-1}e^{-b'λ}} dλ \\
&= \frac{C_G(a',b')}{x_*!} \int{λ^{x_* + a'-1}e^{-(1 + b')λ}} dλ \\
\end{align}

ここで、ガンマ分布の尤度関数
$$ Gam(λ|a, b) = C_G(a, b) λ^{a - 1} e^{- b λ} $$
より、
$$ λ^{a - 1} e^{- b λ} = \frac{Gam(λ|a, b)}{C_G(a, b)}$$
である。

\begin{align}
p(x_*)
&= \frac{C_G(a',b')}{x_*!} \int{λ^{x_* + a'-1}e^{-(1 + b')λ}} dλ \\
&= \frac{C_G(a',b')}{x_*!} \int{\frac{Gam(x_* + a', 1 + b')}{C_G(x_* + a', 1 + b')}} dλ \\
&= \frac{C_G(a',b')}{x_*!} \frac{\int{Gam(x_* + a', 1 + b')} dλ}{C_G(x_* + a', 1 + b')} \\
&= \frac{C_G(a',b')}{C_G(x_* + a', 1 + b') x_*!} \\
\end{align}

ここで、ガンマ分布の尤度関数
$$ Gam(λ|a, b) = C_G(a, b) λ^{a - 1} e^{- b λ} $$
の定数 $ C_G $ を再掲。
$$ C_G(a, b) = \frac{b^a}{Γ(a)} $$
である。

\begin{align}
p(x_*) 
&= \frac{C_G(a',b')}{C_G(x_* + a', 1 + b') x_*!} \\
&= \frac{{b'}^{a'}}{Γ(a')} \frac{Γ(x_* + a')}{(1 + b')^{x_* + a'}} \frac{1}{x_*!} \\
&= \frac{Γ(x_* + a')}{Γ(a')x_*!} (\frac{b'}{1 + b'})^{a'} (\frac{1}{1 + b'})^{x_*} \\
\end{align}

全然ピンとこないが、実は最後の式は、負の二項分布として知られる分布の尤度関数の形をしているらしい。

負の二項分布は、パラメータ $ r \in \mathbb{R^+} $ および $ p \in (0, 1) $ をもつ離散確率分布で、その尤度関数は以下。

\begin{align}
p(x)
&=& NB(x|r, p) \\
&=& \frac{Γ(x + r)}{Γ(r) x!} (1 - p)^r p^x \\
\end{align}

ここで、期待値は $ \frac{pr}{1 - p} $ らしい。

予測分布の計算でよく存じ上げない分布が登場したので、本ページの最後で補足する。

とにかく、共役事前分布を用いると、
ポアソン分布のパラメータ $ λ $ の従うガンマ分布の更新が解析的にできてよかった。

負の二項分布

負の二項分布について。

幾何分布（ベルヌーイ試行が初めて成功するまでの失敗の回数の分布）との関係だったり解釈については省略。

実装と期待値の確認だけ。

\begin{align}
NB(x|r, p) &= \frac{Γ(x + r)}{Γ(r) x!} (1 - p)^r p^x \\
&= \frac{Γ(x + r) p^x}{x!}・\frac{(1 - p)^r}{Γ(r)}
\end{align}

パラメータは $ r \in \mathbb{R^+} $ および $ p \in (0, 1) $

期待値は $ \frac{pr}{1 - p} $

パラメータ $ (r,p) = (2, 0.4), (4, 0.4), (4, 0.6)$ についてグラフを書いてみる。

import numpy as np
from scipy.special import gamma
from scipy.special import factorial

# Negative Binomial Distribution
def NegBin(x, r, p):
  # xに関する項
  main = gamma(x + r) * p ** x / factorial(x)
  # 定数項
  C_NB = (1 - p) ** r / gamma(r)
  return main * C_NB

import matplotlib.pyplot as plt

r_p = [(2, 0.4), (4, 0.4), (4, 0.6)]

# 描画用関数
def show_NegBin(r_p):
  r, p = r_p[0], r_p[1]
  x = np.linspace(start=-10, stop=10, num=2000)
  y = NegBin(x, r, p)
  # (r,p)の組と、期待値('Exp')を表示
  Exp = np.around(r * p /(1 - p), decimals=2)
  plt.title(label=f'(r,p) = {r_p}, Exp = {Exp}', loc='center')
  plt.ylim([0, 1.0])
  plt.plot(x,y)

# 描画するよ
plt.figure(figsize=(12,3))

plt.subplot(131)
show_NegBin(r_p[0])

plt.subplot(132)
show_NegBin(r_p[1])

plt.subplot(133)
show_NegBin(r_p[2])

一旦、ポアソン分布、ガンマ分布、負の二項分布の解釈性については保留しておくが、共役事前分布を用いると解析的に事後分布、そして予測分布が得られることが示せた。

次回は、混合モデルを抽象的に定式化し、ポアソン混合モデルに適用していこうと思う。