はじめに

ベイズ推論による機械学習を読んだので、実装の練習も兼ねていくつかさらっとメモを残しておこうと思います。
今回は、ポアソン混合分布のギブスサンプリングによる推論です。

ポアソン混合分布の推定

ポアソン分布は離散非負の値をとり、以下のヒストグラムは２つのポアソン分布からサンプリングしたものです。 f:id:kenzo1122:20190322234624p:plain

数式では、以下の通りです。

$p(x_n|\bf{s_n}, \bf{\lambda_k})=\prod_k^{K} Poi(x_n|\lambda_k)^{s_{n,k}}$

$\bf{s_n}$ は潜在変数であり、 (0,0,1,0,...,0) のようなクラスタ数Kあるonehot表現されたベクトルです。つまり、 $s_{n,k} =1$ というのは、n番目のデータはk番目の分布からサンプリングされたデータであることを示しています

また、潜在変数はカテゴリ分布に従います。 $\bf{π}$ は混合率と言って、混合分布から観測されたデータに対して、各分布の混合の比率をを示すパラメータです。

$p(\bf{s_n}) = Cat ( \bf{s_n} | \bf{π})$ $where\;\; \sum_k^{K} \pi_k = 1$

そして、ポアソン分布のパラメータ $\bf{ λ}$ に対してポアソン分布の共役事前分布であるガンマ分布を導入しますが、今回はガンマ分布のパラメータ $a,b$ は固定とします。

$p( \lambda_k )=Ga( \lambda_k | a, b)$

ギブスサンプリング

ギブスサンプリングのざっくり流れとしては、適当なパラメータ $x_1^{(0)},x_2^{(1)},x_3^{(1)}$ があるとして、 $x_1^{(1)}$ をサンプルするときすでにサンプルされた値 $x_2^{(1)},x_3^{(1)}$ で分布を条件付けしてサンプリングします。続けて、 $x_2^{(2)}$ を得るために $x_1^{(1)},x_3^{(1)}$ で分布を条件付けして〜というようにサンプリングを繰り返すことで、各パラメータごと更新していく手法です

ポアソン混合分布の生成に関わるパラメータ $\bf{π}, \bf{ λ}$ と潜在変数 $\bf{ S }$ の２つの事後分布に分けてサンプルすることで簡単な確率分布を導出できるのでギブスサンプリングを適用できます。

$\bf{ S } \sim p( \bf{S} | \bf{X}, \bf{ λ} ,\bf{π}) \; , \bf{ λ} , \bf{π} \sim p( \bf{ λ} , \bf{π}| \bf{X}, \bf{S})$

ここからは実装に必要なこれらの確率分布の計算結果だけを載せていきます。

潜在変数の事後分布

$p( \bf{S} | \bf{X},\bf{ λ} ,\bf{π}) \propto p( \bf{X} | \bf{S}, \bf{ λ} )p( \bf{S}| \bf{π}) = \prod_n^{N}p( x_n | \bf{s_n}, \bf{ λ} )p( \bf{s_n}| \bf{π})$

データとパラメータが与えられた元では、条件付き独立が成立し、 $\bf{S}$ を独立にサンプリングが可能になることが分かります。

上記はポアソン分布とカテゴリ分布であり、展開して $s_{n,k}$ でまとめるとカテゴリ分布になります。

$\bf{s_n} \sim Cat( \bf{s_n} | \bf{ \eta_n } )$

$\eta_{n,k} \propto exp( x_n \log \lambda_k - \lambda_k + \log \pi_k )$ $\sum_k^{K} \eta_{n,k} = 1$

パラメータの事後分布

ベイズの定理を用いると、

$p(\bf{ λ} , \bf{π} ,\bf{X}, \bf{S}) \propto p( \bf{X} | \bf{S}, \bf{ λ} )p( \bf{S} | \bf{π} )p( \bf{pi})p( \bf{ λ} )$

ここでも $\bf{X} ,\bf{S}$ が与えられた元では、各パラメータの分布が独立に分解されます。

まず $\lambda_k$ の分布は、上記の $\lambda_k$ に関係ある項に対数を取り、

$\log p( \bf{X} | \bf{S}, \bf{ λ} )p(\bf{ λ} ) = \sum_n^{N} \sum_k^{K} s_{n,k} \log Poi(x_n | \lambda_k) + \sum_k^{K} \log Ga( \lambda_k | a,b)$

展開して、 $\lambda_k$ でまとめてみるとガンマ分布になります。

$\lambda_k \sim Ga( \lambda_k | \hat{a_k}, \hat{b_k})$ $\hat{a_k} = \sum_k^{K}s_{n,k}x_n+a$ $\hat{b_k} = \sum_k^{K}s_{n,k} + b$

$\bf{π}$ の分布は、上記の $\bf{π}$ に関係ある項に対数を取り、

$\log p( \bf{S} | \bf{π})p( \bf{π}) = \sum_n^{N} \sum_k^{K} \log Cat( s_{n,k} | \pi_k ) + \sum_k^{K} \log Dir( \pi_k | \alpha_k)$

展開して、 $\bf{π}$ でまとめてみるとディリクレ分布になります。

$\bf{π} \sim Dir( \bf{π} | \bf{\hat{α}})$ $\hat{\alpha_k} = \sum_n^{N}s_{n,k} + \alpha_k$

実装

学習に使うデータは、 $\lambda_1 = 25 , \lambda_2 = 10$ のパラメータの２つのポアソン分布に従う1次元データ(冒頭のヒストグラム)で、それぞれ150個、250個の合計400個のデータを得ていると仮定します。

x1 = np.random.poisson(lam = 25,size = 150)
x2 = np.random.poisson(lam = 10, size = 250)

ソースコード

以下がコードです.

#ギブスサンプリング

import numpy as np
import matplotlib.pyplot as plt
import math
np.random.seed(1103)

def softmax(X):
    l = []
    for x in X:
        sum_ex = np.sum(np.exp(x))
        l.append([np.exp(i)/sum_ex for i in x])
    return np.array(l)

N = 400

#number of class
K = 2

#paramter of poisson dist
lam = np.ones(K).reshape(1,-1)

#parameter of category dist
pi = np.array([0.5,0.5]).reshape(1,-1)

#latent var
S = np.zeros([N,K])

#parameters of gamma dist
a = np.ones(K)
b = np.ones(K)

#parameters of dirichlet
alpha = np.ones(K)

#train data
x1 = np.random.poisson(lam = 25,size = 150)
x2 = np.random.poisson(lam = 10, size = 250)
X = np.r_[x1,x2].reshape(-1,1)

max_iter = 30

lam_sample,pi_sample =[],[]

for i in range(max_iter):
    
    #S
    ex = X.dot(np.log(lam)) - lam + np.log(pi)
    ex = softmax(ex)

    for n in range(N):
        S[n] = np.random.multinomial(1,ex[n],1)
    
    #λ
    a_hat = X.T.dot(S) + a    
    b_hat = S.sum(axis = 0) + b
    lam = np.random.gamma(a_hat,1/b_hat)
    lam_sample.append(lam[0])
    
    #π
    alpha_hat= S.sum(axis = 0) + alpha
    pi = np.random.dirichlet(alpha_hat)

    pi_sample.append(pi)


fig = plt.figure(figsize=(10,5))
pi_sample = np.array(list((map(list,pi_sample))))
lam_sample = np.array(list((map(list,lam_sample))))

plt.subplot(1,2,1)
for i in range(K):
    pi_m =np.mean(pi_sample[:,i])
    plt.hist(pi_sample[:,i],bins=10,alpha=0.6,label=f"{i}: {round(pi_m,2)}")
plt.xlabel("$π$")
plt.legend()
plt.subplot(1,2,2)
for i in range(K):
    lam_m =np.mean(lam_sample[:,i])
    plt.hist(lam_sample[:,i],bins=10,alpha=0.6,label=f"{i}: {round(lam_m,2)}")
plt.xlabel("$π$")
plt.legend()