ポアソン混合モデルのギブスサンプリング

以下の本の4章を実装していく。

機械学習スタートアップシリーズ ベイズ推論による機械学習入門 (KS情報科学専門書)

まずはギブスサンプリングから。

設定

問題設定としては、観測データの生成過程にポアソン混合モデルを仮定して、 その混合要素である各ポアソン分布のパラメータと混合比の分布が知りたい。

アルゴリズム

ギブスサンプリングのアルゴリズムの導出は、上記教科書に丁寧にかかれているので省略。 結論はP.133 のアルゴリズム4.2にまとめられている。結論のみ書いておく。

実験

コードは末尾にある。

まずは疑似データを生成する。

疑似データは、観測量の数N= 1000とし、がそれぞれ5,20,50の3つのポアソン分布(k=3)から生成されているとする。 混合比は適当に3:6:1とした。

以下の観測量のヒストグラム。

以下が、MAXITER=100で実験したときの iter 20 刻みの状況。青が観測量、赤が予測。

iter=40辺りで多峰性（ここでは3）を捉えている。

まぁそもそも簡単なデータなので、だからどうしたという感じではあるが。

コード

numpy で実装した。

import numpy as np

def fit(X, k, maxiter=100):
    
    N = X.shape[0]
    
    # hyperparams
    lambda_k = np.random.rand(k)
    pi_k = np.random.rand(k)
    a_ini = np.random.rand()
    b_ini = np.random.rand()
    alpha_k_ini = np.random.rand(k)
    
    pred_history = []
    for _iter in range(maxiter):
        # sample s_n from Categorical dist.
        eta_nk = np.exp(np.outer(X, np.log(lambda_k)) - lambda_k + np.log(pi_k))
        z = np.linalg.norm(eta_nk, ord=1, axis=1)
        eta_nk = eta_nk / z[:, np.newaxis]    
        s_n = np.empty((N, k))
        for n in range(N):
            s_n[n] = np.random.multinomial(n=1, pvals=eta_nk[n])

        # sample lambda_k from Gamma dist.
        a_hat = s_n.T@X + a_ini
        b_hat = s_n.sum(axis=0) + b_ini
        lambda_k = np.random.gamma(shape=a_hat, scale=1./b_hat)

        # Sample pi_k from Dirichlet dist.
        alpha_hat = s_n.sum(axis=0) + alpha_k_ini
        pi_k = np.random.dirichlet(alpha_hat)
        
        if _iter == maxiter - 1 or _iter % 20 == 0:
            pred_history.append((lambda_k.copy(), pi_k.copy()))
    
    return pred_history

import matplotlib.pyplot as plt

p0 = np.random.poisson(5, 300) 
p1 = np.random.poisson(20, 600)
p2 = np.random.poisson(50, 100)
X = np.hstack((p0,p1,p2))
plt.hist(X, bins=20)
plt.show()

pred_his = fit(X, k=3)

余談 ガンマ分布の定義

• 須山本と numpy.random.gamma におけるパラメータの表現が少し異なるため、APIを使う時には注意が必要。以下の定義を見比べてみると明らかなように、須山本のa、b は numpy.random.gamma のキーワード引数shape (k)とscale( )にそれぞれ対応している。

  • 須山本の定義は式2.46で与えられていて、次式。

    • 

  • 一方 numpy.random.gamma は 次式

    • 

    • 参考 https://docs.scipy.org/doc/numpy-1.15.0/reference/generated/numpy.random.gamma.html

以上です。

次回は変分推論を実装します。