ポアソン混合モデルの生成モデルからサンプリングする

ベイズモデルにおいては、あるモデルを考えるときに、生成モデル(generative model)を考えます。生成モデルとは、データが生成する過程をモデル化することで、データの分布を表現する手法です。

生成モデルによる定式化を理解していないと、ベイズ手法で問題に合わせたモデルを構築することは難しいため、今回は、あえて基本的な問題であるポアソン混合モデルを通して、ベイズモデルを構築する練習をします。

ポアソン混合モデルの概要

例えば、このようなデータが存在したとします。このデータは、多峰性で、正規分布などの単純な確率分布では表現できなそうですよね。混合分布では、このような分布を、複数のポアソン分布の重ね合わせとして表現します。

例えば、こちらのデータは峰が3つほどあるので、3つのポアソン分布の重ね合わせてできているというように表現します。

例えば、パラメータがそれぞれ、20, 50, 80mのポアソン分布の重ね合わせると次のようになります。実際にはこれらの分布を重ね合わせると次のようになります。

一番最初に提示したデータのヒストグラムの真の分布が、3つのポアソン分布の重ね合わせであるというように仮定したモデルが、混合ポアソンモデルになります。

では、次にこの混合ポアソンモデルを定式化し、実際にサンプリングしてみることにします。

ポアソン混合モデルの定式化

まず、ベイズ的にモデルを構築する場合は、生成モデルを考えます。

まず1つのデータを考えたときに、K個のクラスにどこに所属するかが確率的に決まります。各クラスにk所属する確率は、$\pi_k$で表現します。また、クラスがk個あることを考えるので、$\pi$をkつ並べたk次元ベクトルを$\bm{\pi}$とします。

\begin{equation}
\bm{\pi} = [ \pi_1, \pi_2, \cdots, \pi_k ]
\end{equation}

また、kクラスのうちどれか1つが選ばれる確率を$\pi_k$は表現しているので、$\pi_k$は次のような条件を満たします。

\begin{equation}
\sum_{k=1}^{K} \pi_k = 1
\end{equation}

混合ポアソンモデルでは、この$\bm{\pi}$も確率的に決定されるというように考えます。では、この(1)(2)を満たす確率$\bm{\pi}$をどのように表現するかというと、よくディリクレ分布が使われます。

ディリクレ分布については、こちらの記事をご覧ください。確率分布の形状やパラメータについて解説しています。

ディリクレ分布の定義や性質を可視化して理解する

機械学習の分野を勉強していると、一度は目にすることがあるディリクレ分布 (Dirichlet distribution) カテゴリカル分布や多項分布のパラメータを表現する共役事前分布になる性質もあり、ベイズ推定の分野でも […]

ディリクレ分布を用いることで、$\bm{\pi}$の分布を表現することができます。

\begin{equation}
p(\bm{\pi}) = Dir(\bm{\pi} | \bm{\alpha})
\end{equation}

ここで、ディリクレ分布の形状を決定する$\bm{\alpha}$も確率変数として捉えることができますが、今回は定数として扱います。

ここまでの流れで、1つのデータを考えたとき、あるデータがあるクラス$k$に所属する確率は、$\bm{\pi}$で決定され、$\bm{\pi}$はディリクレ分布によって決定されるという確率モデルまで考えることができました。

ちなみに、$\bm{\pi}$のことを混合比率(mixture ratio)と呼ばれています。

では続いて、$\bm{\pi}$によって、あるデータがクラス$k$に決定されたとします。これを数式で表現してみます。

一般的に、あるn番目のデータを$x_n$と表現し、それがクラス$k$に属することを、混合モデルではよく、$\bm{s_{n}}$で表現します。

$\bm{s_{k}}$は次にような、サイズKのベクトルです。

\begin{equation}
\bm{s_n} = [s_{n1}, s_{n2}, s_{n3}, \cdots, s_{nk}]
\end{equation}

ここで、各要素$s_{nk}$はk番目の要素が選択されているときは1になり、それ以外は0になります。

とはいってもイメージできないと思うので、例えば3番目の要素が選択されているときは$\bm{s_{k}}$は

\begin{equation}
\bm{s_n} = [0, 0, 1 , 0, \cdots, 0, 0]
\end{equation}

と3番目の要素が1になり、それ以外の要素は全て0になります。つまり、これを定式化すると、各要素$s_{nk}$は、

\begin{equation}
\sum_{k=1}^{K}s_{nk} = 1 \\
s_{nk} \in \{ 0, 1 \}
\end{equation}

となります。混合モデルの説明をされるときは、いろんな教科書等では、(6)の記述で解説されていることが多いですが、単純に(5)式のように、選べれている要素が1になり、それ以外の要素が0になるという至極単純なベクトルなので、数式というより、イメージで覚えておくと良いでしょう。

ここで、ポアソン混合モデルにおけるベイズモデルでは、$s_{nk}$も確率的に決定されると考えます。

このような分布はカテゴリカル分布によって、次のように表現できます。

\begin{equation}
p(\bm{s_n} | \bm{\pi}) = Cat(\bm{s_n} | \bm{\pi})
\end{equation}

カテゴリカル分布については、こちらの記事をご覧ください。

カテゴリカル分布について分かりやすく解説

カテゴリカル分布は、ベルヌーイ分布を3次元以上に拡張した分布で、例えるなら、出る目が均等ではない歪んだサイコロを1回投げた時に、どの目が出るか？を表現することができる確立分布です。 ベルヌーイ分布がコインの表と裏のような […]

$\bm{s_{n}}$を行列としてまとめたものをベクトル$\bm{S}$で表現し、これを潜在変数と呼ぶことがあります。

\begin{equation}
\bm{S} = [ \bm{s_1}, \bm{s_2}, \cdots, \bm{s_N}]
\end{equation}

これは、混合モデルから生成されるデータが観測されていても、そのデータがどのモデルから生成されるか潜在的に決まっていても、実際に観測されないということで、この潜在的に決まっていることを強調するためか、よく潜在変数と呼ばれています。

しかし、僕もなぜ潜在変数と呼ばれているのか本当のところはわかりません。ただ、この$\bm{S}$が潜在変数と呼ばれていることくらいは頭の片隅に入れておくと良いでしょう。

ここまで生成モデルを考えるにあたって、潜在変数$\bm{S}$により、あるデータ$x_n$がどのクラス$k$に所属するか、までモデル化することができました。

では$x_n$があるクラス$k$に所属するとしたとき、今回ポアソン混合モデルを仮定しているので、ポアソン分布に従ってデータがランダムに生成されることになります。

ここで、今回はポアソン混合モデルを考えているため、ポアソン分布に限定していますが、これが正規分布であってもガンマ分布であっても、線形モデルであっても構いません。混合モデルの強みは、ある程度モデルに限定されずに、好きにモデルを組み立てられる点が挙げられます。

と、話を戻して、ある$k$番目のクラスのポアソン分布を考えていきます。ポアソン分布の確率密度関数は次のようにかけました。

\begin{equation}
Poi(x | \lambda ) = \frac {\lambda^{x}}{x !} e^{- \lambda}
\end{equation}

(8)のように、ポアソン分布は1つのパラメータ$\lambda$で記述できます。ここでクラス$k$のポアソン分布のパラメータを$\lambda_k$とします。

混合ポアソンモデルでは、この$\lambda_k$が確率的にモデリングします。$\lambda_k$は任意の性の実数値を取るため、このパラメータをガンマ分布で表現することができます。ガンマ分布の解説はこちらの記事をご覧ください。

ガンマ分布とは？期待値や分散などをわかりやすく解説

ガンマ分布(Gamma distribution)は、正の実数 $\lambda \in \mathcal{R}$を確率変数の範囲にもつ、連続確率分布です。 よくベータ分布(Beta distribution)と混同され […]

\begin{equation}
p(\lambda_k) = Gam(\lambda_k | a_k, b_k)
\end{equation}

これで、ガンマ分布から全てのクラス$k$のパラメータ$\lambda_k$得ることができました。

この全ての$\lambda_k$をベクトルで$\bm{\lambda}$とおきます。

\begin{equation}
\bm{\lambda} = [\lambda_1, \lambda_2, \cdots, \lambda_k]
\end{equation}

このようにおくと、(8)で表現した$\bm{S}$と$\bm{\lambda}$を用いると、データの生成モデルは次のように表現でいます。

\begin{equation}
p(x_n | \bm{\lambda} , \bm{S} ) = \prod_{k=1}^{K} Poi(x_n | \lambda_k)^{s_{nk }}
\end{equation}

ここまでで確率的にポアソン分布のパーツを表現することができました。続いて、ポアソン混合分布の生成モデルを考えていきましょう。

混合ポアソンモデルの生成モデル

ここまでの流れを整理すると、混合ポアソンモデルにおけるデータ生成モデルは次にようになりました。

という流れでポアソン混合モデルのデータ生成されることになります。このようにベイズモデルを考えるときは、データの生成モデルを考えることになります。

このポアソン混合分布(possion mixture distribution)のグラフィカルモデル(graphical model)を表現すると次にようになります。

グラフィカルモデルを利用すると、パラメータや確率変数が明確になります。

では、ここまでポアソン混合モデルの生成モデルやグラフィカルモデルを記述してきました。では、実際の確率密度関数は次にように書くことができます。

\begin{equation}
\begin{split}
p(x, \bm{s}, \bm{\lambda}, \bm{\pi}) &= 
p(x | \bm{s}, \bm{\lambda}, \bm{\pi})
p(\bm{s}, \bm{\lambda}, \bm{\pi}) \\
&= 
p(x | \bm{s}, \bm{\lambda}, \bm{\pi})
p(\bm{\lambda}|  \bm{s})
p(\bm{s}|  \bm{\pi}) p(\bm{\pi})
\end{split}
\end{equation}

ポアソン混合モデルからサンプリング

(13)式で、ポアソン混合モデルの生成モデルを記述することができたので、これらの生成モデルからポアソン混合モデルに従うデータをサンプリングしてみます。

まず、今回利用するライブラリをインポートしつつ、シード値を設定しておきます。

import numpy as np
from scipy.stats import dirichlet

np.random.seed(1)

今回は、クラスタ数が$K=3$で、それぞれのポアソン分布のパラメタ$\lambda_1$、$\lambda_2$、$\lambda_3$が次にようなガンマ分布から確率的に選ばれると仮定します。

本来ベイズ推論をする場合は、データからこの$lambda$の分布を決めるパラメータ$a_1, b_1$、$a_2, b_2$、$a_3, b_3$を決定していくのですが、今回は生成モデルを考えているだけなので、この値を仮定してデータを取り出してみることにします。

これらのパラメータから、ランダムに$\bm{\lambda}$を取り出すと次のようになります。

K = 3
a1, b1 = 5, 1
a2, b2 = 4, 2
a3, b3 = 3, 8

lambda1 = np.random.gamma(a1, b1)
lambda2 = np.random.gamma(a2, b2)
lambda3 = np.random.gamma(a3, b3)
lambdas = [lambda1, lambda2, lambda3]

続いて、まず、1つのデータをサンプリングすることを考えていきます。

(11)式における、$p(\bm{\pi})$を決定していきます。これはグラフィカルモデルからも分かるように、これはディリクレ分布からサンプリングできます。

alpha = (1, 1, 1)
pi = np.random.dirichlet(alpha = alpha, size =1)

# => [5.917422657060486, 6.986711710034047, 39.99922182447602]

ディリクレ分布を利用することで、サンプリングすることができました。ちなみに、ディリクレ分布のパラメータ$\bm{\alpha}$も本来のベイズ推論ではデータから学習するパラメータですが、今回は定数をとして扱っていきます。

続いて、この$\bm{\pi}$をカテゴリカル分布に与えることで、$\bm{s_n}$をサンプリングしていきます。

sn = np.random.multinomial(1, pi.reshape(3))
# => array([0, 0, 1])

3番目のクラスタが選ばれたことがわかりました。続いて、このクラスタから$x_n$をサンプリングします。(12)式をそのまま実装しています。

xn = 1
for i, s_nk in enumerate(sn):
    x_nk = np.random.poisson(lambdas[i])
    val *= x_nk ** sn[i]

これで、データを1つサンプリングすることができます。では、これを500点ほど取得してみて、分布がどのようになっているか確かめてみましょう。コードの全文は次のようになっています。

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

np.random.seed(1)

K = 3
a1, b1 = 5, 1
a2, b2 = 4, 2
a3, b3 = 3, 8

lambda1 = np.random.gamma(a1, b1)
lambda2 = np.random.gamma(a2, b2)
lambda3 = np.random.gamma(a3, b3)
lambdas = [lambda1, lambda2, lambda3]

alpha = (1, 1, 1)
pi = np.random.dirichlet(alpha = alpha, size =1)

samplings = []

n_iter = 500
for i in range(n_iter): 
    sn = np.random.multinomial(1, pi.reshape(3))
    xn = 1
    for i, s_nk in enumerate(sn):
        x_nk = np.random.poisson(lambdas[i])
        xn *= x_nk ** sn[i]
    samplings.append(xn)
    
plt.hist(samplings, bins=20)

$N=500$でサンプリングすると、得られた混合ポアソン分布のデータをヒストグラムにするとこのようになりました。

通常のポアソン分布では得られないような、多峰性のあるデータ分布を得ることができました。