ポアソン混合モデル　その4：ポアソン混合モデルによる学習、EMアルゴリズムの適用<実装編-1：データの生成>

今回は、EMアルゴリズムを使ってポアソン混合モデルによる学習を行うことにする。実装編です。

理論編はこちら。

ポアソン混合モデル　その3：ポアソン混合モデルによる学習・推論、EMアルゴリズムの適用<理論編>
今回は、ポアソン混合モデルからデータを生成することにしてみよう。

次回は、この記事で作成したデータに基づいて、ポアソン混合モデルによる推論を行う。

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

ポアソン混合モデルからデータが生成する過程(簡略ver)は以下。

１．パラメータのてきとうな初期化

クラスタの混合割合 $ \textbf{π} $ と $ λ_k　(k = 1, 2, .. , K) $ をてきとうに初期化する。

２．クラスタの割り振り

クラスタ番号を表すベクトル $ \textbf{s}_n (n = 1, 2, .., N) $ が、$ Cat(\textbf{s}_n|\textbf{π}) = \prod_{k=1}^{K} {\pi_k}^{s_{n,k}} $ からサンプルされる。

３．ポアソン分布からのサンプリング

$ x_n (n = 1, 2, .., N) $ が、$ \prod_{k=1}^{K} Poi(x_n|λ_k)^{s_{n,k}} $からサンプルされる ($ s_n $ に対応する一つのポアソン分布からの値が採用される)。

１．パラメータのてきとうな初期化

クラスタの混合割合 $ \textbf{π} $ と $ λ_k　(k = 1, 2, .. , K) $ を適当に初期化する

てきとうに、$ K=3 $ 、 $ π = (\frac{3}{6}, \frac{2}{6}, \frac{1}{6}) $ 、 $ \textbf{λ} = (5, 10, 20) $ としよう。

はい完了。

import numpy as np

K = 3
pi_vec = np.array([3/6, 2/6, 1/6])
lam_vec = np.array([5,10,20])

２．クラスタの割り振り

クラスタ番号を表すベクトル $ \textbf{s}_n (n = 1, 2, .., N) $ が、$ Cat(\textbf{s}_n|\textbf{π}) = \prod_{k=1}^{K} {\pi_k}^{s_{n,k}} $ からサンプルされる。

今回は、お手軽にnumpyの numpy.random.choice メソッドを活用させてもらおう。

N = 100

# (0,1,2)から確率(3/6,2/6,1/6)でN個のデータをサンプリング
S = np.random.choice(a=range(K), size=N, p=pi_vec)

上コードでは、$ s_n \in {0,1,2} $ であり、1 of K なベクトルでなく、クラスタ番号を表すスカラー値としている。

実装上こちらの方が便利そうなのでこのまま進めるが、1 of K なベクトルの作り方を一応紹介しておく。

# クラスタの数と同じ次元の単位行列を作る
I = np.eye(K)

one_hot_S = I[S]

３．ポアソン分布からのサンプリング

$ x_n (n = 1, 2, .., N) $ が、$ \prod_{k=1}^{K} Poi(x_n|λ_k)^{s_{n,k}} $ ($ s_n $ に対応する一つのポアソン分布)からサンプルされる。

これの簡単な実装方法としては、例えば以下があげられる。

・$S_k = \sum_{n=1}^{N} s_{n,k}$　として、各クラスタごとに $ S_K $ 回の試行を行う。

・各 $ s_{n,k} $ に対応するクラスタを選択し、該当するポアソン分布からのサンプリングを実行する。

前者だと秒で済むが、今回は後者を実装してみよう。

from scipy import stats
# 各nについて、3つのポアソン分布のそれぞれからのサンプリングを行っちゃう
mixed_poisson = [stats.poisson.rvs(mu=lam_vec[k], size=N) for k in range(K)]

現時点で、サンプルされたすべてのデータを可視化するとこんな感じ。

plt.hist(mixed_poisson[0], alpha=0.7)
plt.hist(mixed_poisson[1], alpha=0.7)
plt.hist(mixed_poisson[2], alpha=0.7)

ここから、$ s_n $ に基づいて採用する値を絞っていく。

'''気持ちはこれ。内法表記を連チャンで使ってる。
samples = []
for i in range(N):
  for s_n in S:
    data = mixed_poisson[s_n][i]
    samples.append(data)
'''
data = np.array([mixed_poisson[s][i] for i in range(N) for s in S])

まとめて表示するとこんな感じ

plt.hist(data)

まとめて表示してもわかりづらいので、クラスタ別にヒストグラムを描いてみよう。

plt.hist(mixed_poisson[0][np.where(S == 0)], alpha=0.7)
plt.hist(mixed_poisson[1][np.where(S == 1)], alpha=0.7)
plt.hist(mixed_poisson[2][np.where(S == 2)], alpha=0.7)

はい、クラスタによって値の採用率が異なっている感じがわかる。

このデータを使って次の記事ではポアソン混合モデルのパラメータ推定をしたいが、データとして保存したりするのがめんどくさいので、データ生成用コードをまとめておこう。

'''１．パラメータのてきとうな初期化'''
import numpy as np

K = 3
pi_vec = np.array([3/6, 2/6, 1/6])
lam_vec = np.array([5,10,20])

'''２．クラスタの割り振り'''
N = 3000

# (0,1,2)から確率(3/6,2/6,1/6)でN個のデータをサンプリング
S = np.random.choice(a=range(K), size=N, p=pi_vec)

'''３．ポアソン分布からのサンプリング'''
from scipy import stats
# 各nについて、3つのポアソン分布のそれぞれからのサンプリングを行っちゃう
mixed_poisson = [stats.poisson.rvs(mu=lam_vec[k], size=N) for k in range(K)]

data = np.array([mixed_poisson[s][i] for i in range(N) for s in S])

まとめるとあんまり難しい内容もなく簡単な感じ。ライブラリって便利ですね。

今日はここまで。