[ベイズ推論による機械学習入門]多次元ガウス混合モデルにおける推論(ギブスサンプリング編)

須山敦志著の『[ベイズ推論による機械学習入門](https://www.amazon.co.jp/%E6%A9%9F%E6%A2%B0%E5%AD%A6%E7%BF%92%E3%82%B9%E3%82%BF%E3%83%BC%E3%83%88%E3%82%A2%E3%83%83%E3%83%97%E3%82%B7%E3%83%AA%E3%83%BC%E3%82%BA-%E3%83%99%E3%82%A4%E3%82%BA%E6%8E%A8%E8%AB%96%E3%81%AB%E3%82%88%E3%82%8B%E6%A9%9F%E6%A2%B0%E5%AD%A6%E7%BF%92%E5%85%A5%E9%96%80-KS%E6%83%85%E5%A0%B1%E7%A7%91%E5%AD%A6%E5%B0%82%E9%96%80%E6%9B%B8-%E9%A0%88%E5%B1%B1-%E6%95%A6%E5%BF%97/dp/4061538322/ref=sr_1_1?ie=UTF8&qid=1540446806&sr=8-1&keywords=%E3%83%99%E3%82%A4%E3%82%BA%E6%8E%A8%E8%AB%96%E3%81%AB%E3%82%88%E3%82%8B%E6%A9%9F%E6%A2%B0%E5%AD%A6%E7%BF%92)』の第4章を読みました。4.4章 ガウス混合モデルにおける推論 をpythonで実装してみました。全コードは[私のgit](https://github.com/kazucmpt/mixture_model_of_gaussian)にあります。

#記号の使い方

小文字の太字は縦ベクトル、大文字の太字は行列は表します。

###多次元ガウス分布

$D$次元ガウス分布を、平均パラメータベクトル$\boldsymbol{\mu}\in \mathbb{R}^D$, 共分散行列${\bf \Sigma}\in \mathbb{R}^{D \times D}$を用いて
$$
{\mathcal N}({\bf x}|\boldsymbol{\mu},{\bf \Sigma})=\frac{1}{\sqrt{(2\pi)^D|{\bf \Sigma}|}} {\rm exp} \left\\{ -\frac{1}{2}({\bf x}-\boldsymbol{\mu})^{\rm T}{\bf \Sigma^{-1}({\bf x}- \boldsymbol \mu)} \right\\}
$$
で定義します。簡略のために精度行列${\bf \Lambda}={\bf \Sigma}^{-1}$を用います。

###ウィシャード分布

$D\times D$の正定値行列${\bf \Lambda}$を生成するウィシャード分布を$D\times D$の正定値行列${\bf W}$を用いて次で表す。

$$
{\mathcal W}({\bf \Lambda}|\nu,{\bf W})= N |{\bf \Lambda}|^{\frac{\nu-D-1}{2}}{\rm exp}\left( -\frac{1}{2}{\rm Tr}({\bf W}^{-1}{\bf \Lambda}) \right)
$$

$N$は規格化定数で、パラメータ$\nu$は$\nu>D-1$を満たす。

###カテゴリ分布
1 of K ベクトル${\bf s}$を生成するカテゴリ分布を次で表す。
$$
\rm{Cat}({\bf s}|\boldsymbol \pi) = \prod_{k=1}^K \pi_k^{s_k}
$$

ここで$\boldsymbol \pi= (\pi_1,\cdots , \pi_K)^{\rm T}$は分布を決めるパラメータで、$\pi_k \in (0,1)$かつ$\sum_{k=1}^K \pi_k = 1$を満たす。

###ディレクレ分布

ベータ分布を多次元に拡張したディレクレ分布を次で表す。

$$
\rm{Dir}({\boldsymbol \pi}|\boldsymbol \alpha ) = N \prod_{k=1}^K \pi_k^{\alpha_k-1}
$$

ここで$N$は規格化定数で$K$次元のパラメータ$\boldsymbol \alpha$のそれぞれの要素は正の実数値です。ディレクレ分布は$\pi_k \in (0,1)$かつ$\sum_{k=1}^K \pi_k = 1$を満たす$\boldsymbol \pi$を生成する。

#理論
詳細は『[ベイズ推論による機械学習入門](https://www.amazon.co.jp/%E6%A9%9F%E6%A2%B0%E5%AD%A6%E7%BF%92%E3%82%B9%E3%82%BF%E3%83%BC%E3%83%88%E3%82%A2%E3%83%83%E3%83%97%E3%82%B7%E3%83%AA%E3%83%BC%E3%82%BA-%E3%83%99%E3%82%A4%E3%82%BA%E6%8E%A8%E8%AB%96%E3%81%AB%E3%82%88%E3%82%8B%E6%A9%9F%E6%A2%B0%E5%AD%A6%E7%BF%92%E5%85%A5%E9%96%80-KS%E6%83%85%E5%A0%B1%E7%A7%91%E5%AD%A6%E5%B0%82%E9%96%80%E6%9B%B8-%E9%A0%88%E5%B1%B1-%E6%95%A6%E5%BF%97/dp/4061538322/ref=sr_1_1?ie=UTF8&qid=1540446806&sr=8-1&keywords=%E3%83%99%E3%82%A4%E3%82%BA%E6%8E%A8%E8%AB%96%E3%81%AB%E3%82%88%E3%82%8B%E6%A9%9F%E6%A2%B0%E5%AD%A6%E7%BF%92)』の4.4章を読んでください。ここでは、前提と結論のみを書くに留めます。

##前提

$K$個のガウス分布

```math
\begin{align}
p({\bf x_n}|{\bf \mu}_k,{\bf \Lambda}_k) = {\mathcal N}({\bf x_n}|{\bf \mu}_k,{\bf \Lambda}_k^{-1})\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,k=1,2,\cdots,K
\end{align}
```

から$N$個のデータ${\bf X}=\\{{\bf x_1},\cdots,{\bf x}_N\\}$を得る。${\bf x}_n(n=1,2,\cdots,N)$は$D$次元のベクトル。$\boldsymbol{\mu}_k$,${\bf \Lambda}\_k$は$k$番目のガウス分布の平均パラメータベクトルと精度行列であり、$\boldsymbol{\mu}\_k\in \mathbb{R}^D$, ${\bf \Lambda}\_k\in \mathbb{R}^{D \times D}$である。

ここで、各${\bf x}_n$がどのガウス分布によって生成されたのかを表現する$K$次元のベクトル${\bf s}_n$を導入する。${\bf x}\_n$が$k$番目のガウス分布由来であるとき、${\bf s}\_n$の第$k$成分$s\_{n,k}$が$1$で、他の成分は$0$になる。ここで、$\boldsymbol{\mu}=\\{\boldsymbol{\mu}\_1,\cdots,\boldsymbol{\mu}\_K\\}$,$\boldsymbol{\Lambda}=\\{\boldsymbol{\Lambda}\_1,\cdots,\boldsymbol{\Lambda}\_K\\}$とすれば、さらに

```math
\begin{align}
p({\bf x_n}|{\bf s}_n,{\bf \mu}_k,{\bf \Lambda}_k) = \prod_{k=1}^K {\mathcal N}({\bf x_n}|{\bf \mu}_k,{\bf \Lambda}_k^{-1})^{s_{n,k}}
\end{align}
```
以下のアルゴリズムの目的は、${\bf X}$から$\boldsymbol{\mu}$,${\bf \Lambda}$,${\bf S}$を求めることである。ただし${\bf S}=\\{{\bf s}_1,\cdots,{\bf s}\_N\\}$である。

この観測モデルに対応する事前分布として、ガウス分布の共役事前分布であるガウス・ウィシャード分布を採用する。すなわち、ハイパーパラメータ${\bf m} \in \mathbb{R}^D,\beta \in \mathbb{R}^+,{\bf W} \in \mathbb{R}^{D\times D},\nu > D-1$を用いて


```math
\begin{align}
p(\boldsymbol{\mu}_k,\boldsymbol{\Lambda}_k) = {\mathcal N}(\boldsymbol{\mu}_k|\boldsymbol{m},(\beta\boldsymbol{\Lambda}_k)^{-1}){\mathcal W}({\bf \Lambda}_k|\mu,{\bf W})
\end{align}
```

と書ける。なお、それぞれのクラスタの混合比率$\boldsymbol{\pi}=(\pi_1,\cdots,\pi_K)^{\rm T}$より${\bf s}_n$が決定する。


```math
\begin{align}
p({\bf s}_n|\boldsymbol \pi) = {\rm Cat}({\bf s}_n|\boldsymbol \pi)
\end{align}
```

$p(\boldsymbol \pi)$はカテゴリ分布の共役事前分布であるディレクレ分布を用いて

```math
\begin{align}
p(\boldsymbol \pi) = {\rm Dir}(\boldsymbol \pi|\boldsymbol \alpha)
\end{align}
```

と表現される。$\boldsymbol \alpha$はディレクレ分布のハイパーパラメータである。以上の議論をまとめて、モデル全体の設計図は以下の通りである。

```math
\begin{align}
p({\bf X},{\bf S},\boldsymbol \mu,{\bf \Lambda},\boldsymbol \pi) &= p({\bf X}|{\bf S},\boldsymbol \mu,{\bf \Lambda})p({\bf S}|\boldsymbol \pi)p(\boldsymbol{\mu},\boldsymbol{\Lambda})p(\boldsymbol \pi) \\
&= \left\{ \prod_{n=1}^N p({\bf x_n}|{\bf s}_n,{\bf \mu}_k,{\bf \Lambda}_k)p({\bf s}_n|\boldsymbol \pi)\right\} \left\{\prod_{k=1}^K p(\boldsymbol{\mu}_k,\boldsymbol{\Lambda}_k)\right\} p(\boldsymbol \pi)
\end{align}
```

##結論

以下のようなギブスサンプリングアルゴリズムが導出される。

```math
\begin{align}
{\bf s}_n &\sim {\rm Cat}({\bf s}_n|\boldsymbol \eta_n) \\
{\bf \Lambda}_k &\sim {\mathcal W}({\bf \Lambda}_k | \hat{\nu}_k,\hat{{\bf W}}_k) \\
\boldsymbol \mu_k &\sim {\mathcal N}(\boldsymbol \mu_k | \hat{\bf m}_k,(\beta{\bf \Lambda}_k)^{-1}) \\
\boldsymbol \pi &\sim {\rm Dir}(\boldsymbol \pi | \hat{\boldsymbol \alpha}) 

\end{align}
```

ここで

```math
\begin{align}
\eta_{n,k} &\propto {\rm exp}\left(-\frac{1}{2}  ({\bf x}_n-\boldsymbol{\mu}_k)^{\rm T} {\bf \Lambda} ({\bf x}_n-\boldsymbol{\mu}_k) +\frac{1}{2} {\rm log}|{\bf \Lambda}_k| + {\rm log} \pi_k \right) \\
\hat{\alpha_k} &= \sum_{n=1}^N s_{n,k} +\alpha_k \\
\hat{\beta}_k &= \sum_{n=1}^N s_{n,k} + \beta \\
\hat{\nu}_k &= \sum_{n=1}^N s_{n,k} + \nu \\
\hat{{\bf m}}_k &= \frac{\sum_{n=1}^Ns_{n,k}{\bf x}_n+\beta{\bf m}}{\hat{\beta}_k} \\
\hat{{\bf W}}^{-1}_k &= \sum_{n=1}^N s_{n,k} {\bf x}_n{\bf x}_n^{\rm T} + \beta {\bf m} {\bf m}^{\rm T}-\beta_k {\bf m}_k {\bf m}^{\rm T}_k + {\bf W}^{-1}
\end{align}
```

である。

#実装

まずは2次元ガウス混合モデルのギブスサンプリングを実装する。

number_of_cluster$(=K)$種類の2次元ガウス分布からnumber_of_data$(=N)$個のデータを得る関数datasetを用意する。show_initialはTrueでnumber_of_data個のデータを2dplotする。正解となる平均ベクトルexpectation_vectors$(=\\{\boldsymbol \mu_1,\cdots,\boldsymbol \mu_K\\})$と共分散行列convariance_matrixs$(=\\{{\bf \Lambda}^{-1}\_1,\cdots,{\bf \Lambda}^{-1}\_K\\})$を設定する。デフォルトでnumber_of_cluster=3としているが、number_of_clusterを変化させる場合は、適宜、設定する。

```python
def dataset(number_of_data,number_of_cluster,show_initial):

	expectation_vectors = [[] for i in range(number_of_cluster)]

	expectation_vectors[0] = [+12,+6]
	expectation_vectors[1] = [0,-6]
	expectation_vectors[2] = [-13,+7]

	convariance_matrixs = [[[]] for i in range(number_of_cluster)]

	convariance_matrixs[0] = [[3,1],[1,12]]
	convariance_matrixs[1] = [[11,2],[2,9]]
	convariance_matrixs[2] = [[5,0],[0,18]]

	values = [[] for i in range(number_of_cluster)]

	for k in range(number_of_cluster):
		values[k] = np.random.multivariate_normal(expectation_vectors[k],convariance_matrixs[k],int(number_of_data/number_of_cluster))

	data = np.vstack(values)

	if show_initial:
		plt.scatter(data[:,0],data[:,1])
		plt.xlim([-20,20])
		plt.ylim([-20,20])
		plt.xlabel("x",fontsize=16)
		plt.ylabel("y",fontsize=16)
		plt.show()

	return data

```

s_vector(=${\bf s}_n$)をカテゴリ分布からサンプリングする。eta[n][k]は$\eta\_{n,k}$を意味する。

pi_vectorは${\boldsymbol\pi}$、predict_precision_matrixsとpredict_expectation_vectorsは予測された正解${\bf \Lambda},{\boldsymbol \mu}$を意味する。

```python
def sampling_svectors(number_of_data, number_of_cluster, pi_vector, predict_precision_matrixs, predict_expectation_vectors):
	
	eta = np.zeros((number_of_data,number_of_cluster))

	for n in range(0,number_of_data):
		for k in range(0,number_of_cluster):
			eta[n][k] = math.exp(-0.5*(x[n] - predict_expectation_vectors[k] ) @ predict_precision_matrixs[k] @ (x[n] - predict_expectation_vectors[k] ) + 0.5*math.log(np.linalg.det(predict_precision_matrixs[k])) + math.log(pi_vector[k]))

	for n in range(0,number_of_data):
		normalized_const = 0
		for k in range(0,number_of_cluster):
			normalized_const = normalized_const + eta[n][k]

		s_vectors[n] = np.random.multinomial(1,eta[n]/normalized_const)

	return s_vectors
```


次にpi_vector($={\boldsymbol \pi}$)をディレクレ分布からサンプリングする。ここでsum2は$\sum_{n=1}^Ns_{n,k}$を意味する。


```python
def sampling_pi_vector(number_of_data,number_of_cluster,s_vectors):

	alpha_vector = np.ones(number_of_cluster)	
	new_alpha_vector = np.zeros(number_of_cluster)


	for k in range(0,number_of_cluster):

		sum2 = 0

		for n in range(0,number_of_data):
			sum2 = sum2 + s_vectors[n][k]

		new_alpha_vector[k] = sum2 + alpha_vector[k]

	pi_vector = np.random.dirichlet(new_alpha_vector)

	return pi_vector
```

次に、predict_precision_matrixsとpredict_expectation_vectorsをウィシャード分布とガウス分布からサンプリングする。ただし、実際にウィシャード分布からサンプリングするのは大変そうだったので、ウィシャード分布の期待値$<{\bf \Lambda}>=\nu {\bf W}$を採用した。

```python
def sampling_predict_precision_matrixs_and_predict_expectation_vectors(number_of_data,number_of_cluster,s_vectors,):

	#ガウスウィシャード分布の超パラメータ
	beta = 1.0
	v = 2 # > D-1
	m_vector = np.ones(2)
	W_matrix = np.identity(2)

	new_beta = np.zeros(number_of_cluster)
	new_m_vector = np.zeros((number_of_cluster,2))


	for k in range(0,number_of_cluster):
		sum1 = 0
		sum2 = 0

		for n in range(0,number_of_data):
			sum1 = sum1 + s_vectors[n][k]*x[n]
			sum2 = sum2 + s_vectors[n][k]

		new_beta[k] = sum2 + beta 
		new_m_vector[k] = (sum1 + beta*m_vector)/new_beta[k]

	for k in range(0,number_of_cluster):
		sum1 = 0
		sum2 = 0
	
		for n in range(0,number_of_data):
			sum1 = sum1 + s_vectors[n][k]*(x[n].reshape(2,1) @ x[n].reshape(1,2))
			sum2 = sum2 + s_vectors[n][k]
	
		new_W_inverse = sum1 + beta*(m_vector.reshape(2,1) @ m_vector.reshape(1,2)) - new_beta[k]*(new_m_vector[k].reshape(2,1) @ new_m_vector[k].reshape(1,2)) + np.linalg.inv(W_matrix)
		new_v = sum2 + v
	
		predict_precision_matrixs[k] = new_v * np.linalg.inv(new_W_inverse)
		predict_expectation_vectors[k] = np.random.multivariate_normal(new_m_vector[k],np.linalg.inv(new_beta[k]*predict_precision_matrixs[k]))

	return predict_precision_matrixs,predict_expectation_vectors
```

二次元ガウス分布を返す変数を用意しておく。

```python
def multivariate_Gaussian_distribution(x,y,precision_matrix,expectation_vector):

	cordinate = np.array([x,y])
	convariance_matrix = np.linalg.inv(precision_matrix)
	
	return 1.0/((2.0*math.pi)**2*np.linalg.det(convariance_matrix))*math.exp(-0.5*(cordinate-expectation_vector)@np.linalg.inv(convariance_matrix)@(cordinate-expectation_vector))

```

求めた$\bf s$を用いてクラスタリングし、散布図をカラーリングして、等高線をひく関数を用意する。

```python
def show_answer(number_of_data,number_of_cluster,x,s_vectors,predict_precision_matrixs,predict_expectation_vectors,show_contour):
	
	groups = [[] for i in range(number_of_cluster)] #空の配列をnumber_of_cluster個用意する。

	for n in range(0,number_of_data):
		for k in range(0,number_of_cluster):
			if s_vectors[n][k] == 1:
				groups[k].append(x[n].tolist())

	for k in range(0,number_of_cluster):
		groups[k] = np.array(groups[k])
		plt.scatter(groups[k][:,0],groups[k][:,1])
	
	if show_contour:

		N = 60
		cordinate_x = np.linspace(-20,20,N)
		cordinate_y = np.linspace(-20,20,N)
		cordinate_X,cordinate_Y = meshgrid(cordinate_x,cordinate_y)

		value_of_z = np.zeros((N,N))

		for k in range(0,number_of_cluster):
			for i in range(0,N):
				for j in range(0,N):
					value_of_z[j][i] = multivariate_Gaussian_distribution(cordinate_x[i],cordinate_y[j],predict_precision_matrixs[k],predict_expectation_vectors[k])

			plt.contour(cordinate_X,cordinate_Y,value_of_z) 
	
	plt.ylim([-20,20])
	plt.ylim([-20,20])
	plt.xlabel("x",fontsize=16)
	plt.ylabel("y",fontsize=16)
	plt.show()

```

最後にmain関数を用意する。


```python
if __name__  == "__main__":
	number_of_cluster = 3 #クラスタの数
	trial_times = 80 #サンプリングを行う回数
	number_of_data = 600 #入力データ数

	assert number_of_data % number_of_cluster == 0, "入力データ数をクラス数の倍数にしてください。"

	show_initial = True #分類前の２次元データ分布を見るか否か。
	show_contour = True #分類後に正解等高線を表示するか否か。

	#初期化
	predict_expectation_vectors = np.zeros((number_of_cluster,2))
	predict_precision_matrixs = np.zeros((number_of_cluster,2,2))
	for i in range(0,number_of_cluster):
		predict_precision_matrixs[i] = np.identity(2)

	pi_vector = np.ones(number_of_cluster)
	s_vectors = np.zeros((number_of_data,number_of_cluster))

	#データを渡す
	x = dataset(number_of_data,number_of_cluster,show_initial)

	#ギブスサンプリングを実行
	for i in range(0,trial_times):

		#sをサンプリングする。
		s_vectors = sampling_svectors(number_of_data, number_of_cluster, pi_vector, predict_precision_matrixs, predict_expectation_vectors)
		
		#predict_precision_matrixsとpredict_expectation_vectorsをサンプリングする。
		predict_precision_matrixs, predict_expectation_vectors = sampling_predict_precision_matrixs_and_predict_expectation_vectors(number_of_data, number_of_cluster, s_vectors)

		#pi_vectorをサンプリングする。
		pi_vector = sampling_pi_vector(number_of_data, number_of_cluster, s_vectors)

	#答えをプロット
	show_answer(number_of_data,number_of_cluster,x,s_vectors,predict_precision_matrixs,predict_expectation_vectors,show_contour)

	#最終結果を表示
	print("Trial回数:",trial_times,"\n")
	for k in range(0,number_of_cluster):
		print("グループ"+str(k)+"の平均\n",predict_expectation_vectors[k])
		print("グループ"+str(k)+"の共分散行列\n",np.linalg.inv(predict_precision_matrixs[k]),"\n")
```

#結果

##2次元ガウス混合モデル

1枚目の画像の600点を3つのクラスタに分類した後の様子が2枚目である。80回で正しくクラスタリングされている様子が分かる。

![Figure_1.png](https://qiita-image-store.s3.amazonaws.com/0/139546/e2105e42-bd12-0ef0-7961-a7d9f523e21c.png)
![Figure_1-1.png](https://qiita-image-store.s3.amazonaws.com/0/139546/15565b03-7c5b-f672-ba05-537ab2672ae4.png)


>グループ0の平均
 [-13.03527306   6.85378566]
グループ0の共分散行列
 [[ 5.55089831 -0.06026977]
 [-0.06026977 20.82277496]] 

>グループ1の平均
 [-0.46999732 -6.09191948]
グループ1の共分散行列
 [[9.62733973 1.28850163]
 [1.28850163 9.60439715]] 

>グループ2の平均
 [12.0153539   5.74785767]
グループ2の共分散行列
 [[3.78214067 1.07161246]
 [1.07161246 9.74919818]] 

##3次元ガウス混合モデル
3次元版でもやってみた！（ソースコードは[こちら](https://github.com/kazucmpt/mixture_model_of_gaussian/blob/master/mixture_model_of_gaussian3D.py)）

![Figure_1-3.png](https://qiita-image-store.s3.amazonaws.com/0/139546/2da961d2-2084-c248-1f23-e9568f201bbc.png)
![Figure_1-2.png](https://qiita-image-store.s3.amazonaws.com/0/139546/34d56af1-9ac0-abe8-aadc-2431f99c50c4.png)