EMアルゴリズムのサンプル

データの分類

「これなら分かる最適化数学」p145
について、次のような問題ありました。
一次元なので次は多次元でやります。ええ。

[改・例題5.4]
それぞれ平均 $\mu_1$ 、 $\mu_2$ 分散 $\sigma_1^2$ 、 $\sigma_2^2$
の正規分布から独立に発生した $N_1,N_2$ 個のデータが総数 $N=N_1+N_2$ 個のデータ{ $x_{\alpha}$ }, $\alpha = 1, \dots , N$ を観測したとする。
これをEMアルゴリズムでクラスを判別せよ。

EMアルゴリズムは以下のようなものです。

所属確率 $w_{\alpha}^{(k)},k=1,\dots,K,\alpha = 1, \dots, N$ の初期値を与える。

次の計算を{ $w_{\alpha}^{(k)}$ }が収束するまで反復する
1.　次のように $N_k,k=1,\dots,K$ を計算する。
$N_k = \sum_{a=1}^{N} w_{\alpha}^{(k)}$

2.次のように $\mu_{\alpha},\sigma_k^2,k=1,\dots,K$ を計算する。
$\mu_k =\frac{1}{N_k}\sum_{a=1}^{N} w_{\alpha}^{(k)}x_{\alpha}$
$\sigma_k^2 =\frac{1}{N_k}\sum_{a=1}^{N} w_{\alpha}^{(k)}(x_{\alpha}-\mu_k)^2$

3.次のように確率密度 $p_k(x),k=1,\dots,N$ を定義する。
$p_k(x)=\frac{1}{\sqrt{2 \pi \sigma_k^2}} exp(-\frac{(x-\mu_k)^2}{2 \sigma_k^2})$

4.次のように所属確率 $w_{\alpha}^{(k)},k=1,\dots,K,\alpha=1,\dots,N$ を更新する。
$w_{\alpha}^{(k)}=\frac{\pi_k p_k(x_{\alpha})}{\sum_{l=1}^2\pi_l p_l(x_{\alpha})}$

$\pi_k = \frac{N_k}{N}$

各 $x_{\alpha}$ を $w_{\alpha}^{(k)},k=1,\dots,K$ が最大になるクラスkに分類する。

そして、以下がそのソースコード。
今回は小学生と成人を混ぜたデータで分類させた。
また初期値は、それぞれのデータの値がすべてのデータの中の大きい順に並べたときの
順位を使った。例えば、もっとも高い数値をもつデータはデータの総数のうち1/100となり、
0.99の確率で大人、0.01の確率で小学生として初期化してある。

また誤り率が出力されるが、20%以下だった。
まぁデータによってかわるんだろうけど。

#-*- coding:utf-8 -*-
import numpy as np
from scipy import stats as sp

Trial_Number=1000 
N1 = 70.
N2 = 30.
N = N1+N2

Height_1 = 170.
Height_2 = 145.

variance_1 = 5.
variance_2 = 10.

#所属確率の初期値を与える
def Initial_W(xline):
	zline = []	
	for i in xline:
		for j in range(len(xline)):
			if i == np.sort(xline)[j]:
				zline.append( [len(xline) - j, j ])
	return np.array(zline).T


#アルゴリズムの2
def E_step(Z,X):
	N_1 = sum(Z[0])
	N_2 = sum(Z[1])
	U_1 = sum(Z[0]*X)/N_1
	U_2 = sum(Z[1]*X)/N_2
	variance_1 = sum( Z[0]*(X-U_1)**2 )/N_1
	variance_2 = sum( Z[1]*(X-U_2)**2 )/N_2
	return N_1,N_2,U_1,U_2,variance_1,variance_2


#アルゴリズムの3(確率密度)
def P_k(xline,expectation,variance):
	return 1/np.sqrt(2*np.pi*variance) *(np.e)**( -(xline-expectation)**2/( 2*variance ) )

#所属確率の更新
def W_update(wline,xline):
	N_1,N_2,U_1,U_2,variance_1,variance_2=E_step(wline,xline)
	PI_1 = N_1/N
	PI_2 = N_2/N
	for i in range(len(xline)):
		P_1 = P_k(xline[i],U_1,variance_1)
		P_2 = P_k(xline[i],U_2,variance_2)
		wline[0][i] = PI_1*P_1/np.sum(PI_1*P_1+PI_2*P_2)
		wline[1][i]=1-wline[0][i]
	return wline


if __name__=="__main__":
	#adults
	A=[ sp.norm.rvs(Height_1,variance_1) for i in range(N1) ]
	#elementary
	ES=[ sp.norm.rvs(Height_2,variance_2) for i in range(N2) ]

	#bind data
	total = A+ES
	
	zline=Initial_W(total)/N

	for i in range( Trial_Number ):
		W_update(zline,total)

	print 1*(zline[0]>.5)
	print "誤り率",(sum(zline[0][:N1] >.5)+sum(zline[1][N1:] > 0.5))/N