概要
ベイズ推論による機械学習入門の内容について、実際にpythonで試してみた際の記録です。
ポアソン混合モデルをガウスサンプリングを使って推論する箇所について記載しました。
実施内容
以下のように作成したデータから、ポアソン混合モデルの事後分布を推論します。
import numpy as np
import matplotlib.pyplot as plt
np.random.normal(0, 1, 10)
poi_test = np.random.poisson([10, 20],(1000,2))
s_test = np.random.multinomial(1, [0.5,0.5], (1000))
x = np.prod(np.power(poi_test, s_test), axis=1)
plt.hist(x,bins=40, range=(0,40))
ポアソン混合モデルから出力したデータから、ポアソン混合モデルの事後分布を推論するのであまり面白くはないですけどねw
ポアソン混合モデルのガウスサンプリング
ポアソン混合モデルのガウスサンプリングは以下のようになります。
sのサンプリング
s_n \sim Cat(s_n|\eta_n) \\
\eta_{n,k} \propto exp\{x_n\ln\lambda_k-\lambda_k + \ln\pi_k\} \\
\left( s.t. \sum_{k=1}^K\eta_{n,k}=1\right)
# sのサンプリング
def sample_s(x, lam, pi):
s = []
for n in range(len(x)):
eta_n = np.exp( x[n] * np.log(lam) - lam + np.log(pi))
eta_n = eta_n / np.sum(eta_n)
s_n = np.random.multinomial(1, eta_n)
s.append(s_n)
s = np.array(s)
return s
λのサンプリング
\lambda \sim Gam(\lambda_k|\hat{a}_k, \hat{b}_k) \\
\hat{a}_k = \sum_{n=1}^Ns_{n,k}x_n + a \\
\hat{b}_k = \sum_{n=1}^Ns_{n,k} + b
# lambdaをサンプリング
def sample_lam(x, s):
a_hat = []
b_hat = []
for i in range(2):
a_hat_k = np.sum(s[:, i] * x, axis=0) + a
b_hat_k = np.sum(s[:, i], axis=0)
a_hat.append(a_hat_k)
b_hat.append(b_hat_k)
lam = np.random.gamma(np.array(a_hat), 1/np.array(b_hat))
return lam
πのサンプリング
\pi \sim Dir(\pi|\hat{\alpha}) \\
\hat{\alpha}_k = \sum_{n=1}^Ns_{n,k} + \alpha_k
# piをサンプリング
def sample_pi(s):
alpha_hat = np.sum(s, axis=0) + alpha
pi = np.random.dirichlet(alpha_hat)
return pi
サンプリングの実行
# 初期値
a = 2
b = 2
alpha = np.ones((2))
lam = np.random.gamma(a,b,2)
pi = np.random.dirichlet(alpha)
s_samples=[]
lam_samples=[]
pi_samples=[]
for i in range(5):
print(i)
for j in range(100):
s = sample_s(x, lam, pi)
s_samples.append(s)
lam = sample_lam(x, s, a, b)
lam_samples.append(lam)
pi = sample_pi(s, alpha)
pi_samples.append(pi)
サンプルした値をみる
本にはサンプル方法は書かれていたのですが、サンプルしたものをどうすればいいのかは良くわかりませんでした(理解度の問題かもしれませんが)。
そこで色々グラフ化してみます。
サンプルされた値からクラスを決めてみる
だいぶ雑ですが、クラス毎にサンプルされた値を全て足して、大きい方にクラス分けします。こういう方法で分けるのが正しいのかはわかりません。。。
x_1 = []
x_2 = []
tmp = np.sum(np.array(s_samples),axis=0)
for i in range(len(x)):
if tmp[i][0] > tmp[i][1] :
x_1.append(x[i])
else:
x_2.append(x[i])
10回サンプルした場合
10回くらいだとあまりいい感じに別れませんね。