k-meansアルゴリズム書いた。

PRMLの定義に従ってk-meansアルゴリズム書いたった。
汚いのであとで整形する予定
とりあえずn次元のクラスタに分けることが可能。

import numpy as np
import pylab as plt
#クラスタの数
K=2 

data=np.genfromtxt("faithful.txt")
X=standardizing(data)

#分けたいクラスターの数K
def k_means(X,k):
	CN=100
	mu,r=initialize(X,k)
	for i in range(CN):
		old_r=r.copy()
		E_step(X,mu,r)
		if np.all(r==old_r): break
		M_step(X,mu,r)
		print i
	return r
	
#initializing mu,r
def initialize(X,k):
	N=data.shape[0]
	K=data.shape[1]
	mu=np.random.rand(k,K)
	r=np.zeros((N,k))
	return mu,r


#n番目のデータ点を中心がそれに最も近いクラスター中心に割り当てる
def E_step(X,mu,r):
	N=r.shape[0]
	K=r.shape[1]
	for n in range(N):
		for k in range(K):
			if k==np.argmin(((X[n]-mu)**2).sum(axis=1)):
				r[n][k]=1
			else:
				r[n][k]=0
	return r

#標準化
def standardizing(data):
	mu=data.mean(axis=0)
	sigma=data.std(axis=0)
	return	(data-mu)/sigma


#目的関数の定義
def J(X,mu,r):
	N=r.shape[0]
	K=r.shape[1]
	J=0.0	
	for n in range(N):
		for j in range(K):
			J+=r[n, j]*np.linalg.norm(X[n]-mu[j])**2
 	return J

#M_step
def M_step(X,mean,r):
	N=mean.shape[0]
	K=mean.shape[1]
	A=np.dot(r.T,X)
	B=r.sum(axis=0)
	for n in range(N):
		for k in range(K):
			mean[n][k]=A[n][k]/B[n]
	return mean

argmin() 関数とても便利.....
"Numpy Example List"

in [78]: import numpy as np
In [79]: a = np.array([10,20,30])

In [80]: a
Out[80]: array([10, 20, 30])

In [81]: minindex = a.argmin()

In [82]: minindex
Out[82]: 0

In [83]: a[minindex]
Out[83]: 10