Pythonで混合モデルの教師なし学習アルゴリズムを実装しました。

www.amazon.co.jp

 この本の第5章(統計的最適化)のクラス分類問題のところです。
 \{x_\alpha\},\alpha=1,...,Nというデータを、教師データを使わずに、k個のクラス(ここでは2クラス)に分類します。各データx_\alphaが両方のクラスにまたがってある確率で所属していると考え、w_\alpha^{(k)}はデータx_\alphaがクラスkに属する確率とします。初期値の与え方に工夫が必要です。

import numpy as np
import random

def unsupervised(N=20,n=10,iteration=1000): #N is number of samples
	#normal distribution
	x1 = np.random.normal(10,5,n)
	x2 = np.random.normal(20,10,N-n)
	x = np.append(x1,x2)
	#initial value of w, vector length is N	
	w1 = np.ones(N)*0.2
	w1[x<15] = 0.8 # for class 1
	w2 = 1 - w1  # for class 2

	for i in range(iteration):
		#calculate N, mu, sigma for each k = 1, 2
		N1 = sum(w1)
		N2 = sum(w2)
		mu1 = sum(w1*x) / N1
		mu2 = sum(w2*x) / N2
		sigma1sq = sum(w1*(x - mu1)**2) / N1
		sigma2sq = sum(w2*(x - mu2)**2) / N2
		#calculate pi and p
		pi1 = N1 / N
		pi2 = N2 / N
		p1 = (1/np.sqrt(2 * np.pi * sigma1sq)) * np.e ** ( -(x - mu1)**2 / (2*sigma1sq) )
		p2 = (1/np.sqrt(2 * np.pi * sigma2sq)) * np.e ** ( -(x - mu2)**2 / (2*sigma2sq) )
		#calculate the new weight
		w1_new = pi1 * p1 / (pi1 * p1 + pi2 * p2)
		w2_new = pi2 * p2 / (pi1 * p1 + pi2 * p2)
		#convergence
		if sum ( abs(w1_new - w1) + abs(w2_new - w2) ) < 0.000001:
			break
		w1 = w1_new
		w2 = w2_new
	return w1, w2, i