ポアソン混合モデルを使ったクラスタリング

# 概要
[ベイズ推論による機械学習入門](https://www.kspub.co.jp/book/detail/1538320.html)の内容について、実際にpythonで試してみた際の記録です。

ポアソン混合モデルをガウスサンプリングを使って推論する箇所について記載しました。

# 実施内容
以下のように作成したデータから、ポアソン混合モデルの事後分布を推論します。

```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))
```

![sample_data2.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/223111/8c34f4a0-1f92-fc01-d6de-ae2f13d97677.png)


ポアソン混合モデルから出力したデータから、ポアソン混合モデルの事後分布を推論するのであまり面白くはないですけどねw

# ポアソン混合モデルのガウスサンプリング
ポアソン混合モデルのガウスサンプリングは以下のようになります。

## sのサンプリング

```math
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)
```

```python
# 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
```

## λのサンプリング

```math
\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
```

```python
# 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
```

## πのサンプリング

```math
\pi \sim Dir(\pi|\hat{\alpha}) \\
\hat{\alpha}_k = \sum_{n=1}^Ns_{n,k} + \alpha_k
```

```python
# piをサンプリング
def sample_pi(s):
  alpha_hat = np.sum(s, axis=0) + alpha
  pi = np.random.dirichlet(alpha_hat)
  return pi
```

# サンプリングの実行

```python
# 初期値
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)
```

# サンプルした値をみる
本にはサンプル方法は書かれていたのですが、サンプルしたものをどうすればいいのかは良くわかりませんでした（理解度の問題かもしれませんが）。
そこで色々グラフ化してみます。

## サンプルされた値からクラスを決めてみる
だいぶ雑ですが、クラス毎にサンプルされた値を全て足して、大きい方にクラス分けします。こういう方法で分けるのが正しいのかはわかりません。。。

```python
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])
```

正解のデータはこんな感じ
![origin_class.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/223111/c907f3fd-1ab9-ec59-42be-2c5d5bc3a74e.png)


### 10回サンプルした場合
10回くらいだとあまりいい感じに別れませんね。

![sample_10.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/223111/43159181-8428-81d7-5d2a-1b91a46a4106.png)


## 20回サンプルした場合
![sample_20.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/223111/81a5052a-49fb-2c28-ad7d-cec82311a68a.png)

## 1000回サンプルした場合
なんだか思いの外バツっと別れました。
![sample_1000.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/223111/dffc99d2-37f9-0780-9783-f8fd81e191de.png)


## lambdaの分布をみてみる
これは最初に決めた10,20あたりにピークが来ました。ちゃんと推論できたってことでいいのかな？
![lambda.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/223111/7d7ad767-5bbb-6fcd-8798-3256a3e2cbb0.png)