今回は前回、理論を確認したK-means法の実装をPythonで行います。
www.wantanblog.com
データ準備
まずはクラスタリングを行うデータを準備します。
対象データを他に準備している場合には、参考程度に無視してぶっ飛ばしてOK。
データセット
データは毎度おなじみ「scikit-learn」のアヤメの品種データ（Iris plants dataset）を扱います。
以下、今回使用する範囲で簡単に振り返りますが、詳細については過去記事もご参照ください。
■説明変数
| 0 | sepal length (cm) | がく片の長さ |
|---|---|---|
| 1 | sepal width (cm) | がく片の幅 |
■目的変数
| 0 | setosa（ヒオウギアヤメ） |
|---|---|
| 1 | versicolor（ブルーフラッグ） |
| 2 | virginica（バージニカ） |
上記のデータを扱い、「がく片の長さ」（X軸）と「がく片の幅」（Y軸）から三種の分類を行ってみます。
今回は「教師なし学習」を扱っているため本来は答えのないデータを扱うべきですが、データセットには答えをもっているためある意味検証はやりやすいです。
データの取得
データをダウンロードしファイルに格納します。
# データ準備 import numpy as npy import matplotlib.pyplot as plt %matplotlib inline # アヤメ品種データの読み込み from sklearn.datasets import load_iris iris_data = load_iris() # 説明変数 X_array = iris_data.data # カラムデータを取得 X0=X_array[:,0] X1=X_array[:,3] # 目的変数 t_array = iris_data.target T=t_array N=len(X0) petal_data=npy.zeros((N,2)) T3 = npy.zeros((N,3), dtype=npy.uint8) #データの設定 for i in range(N): petal_data[i]=[X0[i], X1[i]] if T[i] == 0: T3[i]=[1,0,0] elif T[i] == 1: T3[i]=[0,1,0] else: T3[i]=[0,0,1] X_range0=[min(X0)*0.9,max(X0)*1.1] X_range1=[min(X1)*0.9,max(X1)*1.1] # データをclassdata3.npzファイルに保存する npy.savez('neural_rawdata.npz',X=petal_data,T3=T3,X_range0=X_range0,X_range1=X_range1,X_n=N)
データの準備
・データの分割
ここでは合計150個のデータを半分に分けて訓練データとテストデータに分けています。
必須ではありません。
# 生データファイルから取り出す sample_data = npy.load('neural_rawdata.npz') # 入力値の設定 X=sample_data['X'] # がく片の長さの表示範囲設定 X_range0=sample_data['X_range0'] # がく片の幅の表示範囲設定 X_range1=sample_data['X_range1'] # クラス（答え）の設定 T3=sample_data['T3'] # 訓練データとテストデータに分割する X_test=npy.r_[X[:25,:],X[50:75,:],X[100:125,:]] X_train=npy.r_[X[25:50,:],X[75:100,:],X[125:150,:]] T_test=npy.r_[T3[:25,:],T3[50:75,:],T3[100:125,:]] T_train=npy.r_[T3[25:50,:],T3[75:100,:],T3[125:150,:]] X_col = ['cornflowerblue','black','white'] npy.savez('unsu_data.npz',X_train=X_train, T_train=T_train,X_test=X_test,T_test=T_test,X_range0=X_range0,X_range1=X_range1,X_col=X_col)
・表示範囲の設定
import math # ファイルデータから取り出す file_data = npy.load('unsu_data.npz') # 入力値の設定 X=file_data['X_train'] # がく片の長さの表示範囲設定 X_range0=file_data['X_range0'] X_range0[0]=math.floor(X_range0[0]) X_range0[1]=math.ceil(X_range0[1]) # がく片の幅の表示範囲設定 X_range1=file_data['X_range1'] X_range1[0]=math.floor(X_range1[0]) X_range1[1]=math.ceil(X_range1[1]) # カラムの色情報設定 X_col = file_data['X_col'] N=len(X) K=3
・データサンプルの図示
教師なし学習には必要ありませんが、データが正しく分割できているか図示して確認しておきます。
import matplotlib.pyplot as plt %matplotlib inline # データの表示 def Show_data(x,t): wk,n=t.shape for i in range(n): plt.plot(x[t[:,i]==1,0],x[t[:,i]==1,1],linestyle='none',marker='o',markeredgecolor='black',color=X_col[i],alpha=0.8) plt.grid(True) plt.figure(1,figsize=(8,3.7)) plt.subplot(1,2,1) Show_data(X_train,T_train) plt.xlim(X_range0) plt.ylim(X_range1) plt.title('Training Data') plt.subplot(1,2,2) Show_data(X_test,T_test) plt.xlim(X_range0) plt.ylim(X_range1) plt.title('Test Data') plt.show()
出力結果
K-means法の実装
ここからが本題のK-means法に関する実装になります。
初期値の設定
まず最初に中心ベクトルとクラス指示変数の初期値を設定してしまいます。
# 変数の準備と初期化 # 中心ベクトルの初期値 Mu = npy.array([[5,0.5],[5,1],[5,2]]) # クラス指示変数の初期値 R = npy.c_[npy.ones((N,1),dtype=int),npy.zeros((N,2),dtype=int)]
3つのクラスターの中心ベクトル（Mu）は理論上はどの位置に設定してもよいはずですが、プロットを散布状況を見て上記のように設定。
今回の場合で言えば、最終的なクラスターの位置は概ね予想がつきそうなものですが、予め予想した位置に近い方がいいのか、わざと離した方がいいのかは分かりませんね。
クラス指示変数（R）は全てクラスター０に所属させる形としました。
■グラフ上に図示
def show_prm(x,r,mu,col): for k in range(K): # データ分布の表示 plt.plot(x[r[:,k] == 1,0], x[r[:,k] == 1,1],marker='o',markerfacecolor=X_col[k], markeredgecolor='k',markersize=6,alpha=0.5,linestyle='none') plt.plot(mu[k,0],mu[k,1],marker='*',markerfacecolor=X_col[k], markersize=15,markeredgecolor='k',markeredgewidth=1) plt.xlim(X_range0) plt.ylim(X_range1) plt.grid(True) plt.figure(figsize=(4,4)) R = npy.c_[npy.ones((N,1)),npy.zeros((N,2))] show_prm(X,R,Mu,X_col) plt.title('Mu and R') plt.show()
データプロットと初期値を図示しています。
クラス０は青、クラス１は黒、クラス２は白に設定しているので、ここでは全てのデータは青で表示されています。
クラス指示変数Rの更新
クラス指示変数Rの更新処理を実装します。
# rの更新 def step1_kmeans(x0,x1,mu): N = len(x0) r = npy.zeros((N,K)) for n in range(N): wk = npy.zeros(K) for k in range(K): # 各プロットと各中心ベクトルとの距離を算出 wk[k] = (x0[n]-mu[k,0])**2 + (x1[n]-mu[k,1])**2 # 最小値となるクラスターに所属させる r[n,npy.argmin(wk)]=1 return r plt.figure(figsize=(4,4)) R = step1_kmeans(X[:,0],X[:,1],Mu) show_prm(X,R,Mu,X_col) plt.title('step1') plt.show() print(R)
step1_kmeans()関数は、相互的に更新を行うステップの１つになるのでこの後も使い回します。
Rの更新処理では、 各プロットと各中心ベクトルとの距離を算出し、それを元に最小値となったクラスターに更新します。
npy.argminでは、渡された配列の最小値を判別しそのインデックスを返却します。
処理の確認としての出力結果は以下の通りになりました。
・出力結果
この時点で一番近いクラスターの中心ベクトルに対してクラスリングされている様子が分かると思います。
クラスターの中心ベクトルμの更新
中心ベクトルμの更新処理を実装します。
# 中心ベクトルμの更新 def step2_kmeans(x0,x1,r): mu=npy.zeros((K,2)) for k in range(K): # x0とx1の平均を算出 mu[k,0]=npy.sum(r[:,k]*x0)/npy.sum(r[:,k]) mu[k,1]=npy.sum(r[:,k]*x1)/npy.sum(r[:,k]) return mu # 処理 plt.figure(figsize=(4,4)) mu=step2_kmeans(X[:,0],X[:,1],R) show_prm(X,R,mu,X_col) plt.title('Step2') plt.show()
中心ベクトルμの更新は現在そのクラスターに所属しているデータプロットのx0とx1の平均を算出することによって、重心を求め更新します。
これを全体の中のステップ２とします。
・出力結果例
今回はデータがかなり単純なので一回の更新でかなりそれっぽい位置に配置されています。
μとRの更新処理を任意の回数だけ繰り返す
K-means法はクラス指示関数Rとクラスターの中心ベクトルμを繰り返す手法なので、前述の処理をうまく繰り返す実装を行います。
plt.figure(1,figsize=(10,6.5)) max_it=6 for it in range(0,max_it): plt.subplot(2,3,it+1) R=step1_kmeans(X[:,0],X[:,1],Mu) show_prm(X,R,Mu,X_col) plt.title("{0:d}".format(it+1)) plt.xticks(range(int(X_range0[0]),int(X_range0[1])), "") plt.yticks(range(int(X_range1[0]),int(X_range1[1])), "") Mu=step2_kmeans(X[:,0],X[:,1],R) plt.show()
具体的な処理は前述部分で実装されていますので、今回はそれらの処理を6回繰り返しを行っています。
・出力結果例
最初にデータをプロットしただけの状態と比較してもうまくクラスタリングが行われているように見えます。
プラスα
クラスタリングの結果を検証する
「教師なし学習」ではそもそもクラスタリングの正解が用意されていないケースに使用するため、本来はクラスタリングがうまくできているかということは単純に検証できることではないとは思います。
今回は学習のために「教師あり学習」でも使用した言わば答えがある状態のデータを扱っているのでここで簡単に検証まで行います。
・実装例
# 正解データ=T3 #クラスタリング結果=R anslist=npy.r_[T3[25:50,:],T3[75:100,:],T3[125:150,:]] i = 0 correct_rate = [] for answer in anslist: correct_rate.append((R[i] == answer).all()) i = i + 1 print(correct_rate.count(True)/len(correct_rate))
・出力結果
0.84
結果として84％はうまくクラスタリングできたという結果になりました。
今回のような単純なデータにしては正答率（教師なし学習において正しい言葉か分かりませんが）が低いという印象を受けました。
これはクラスタリングを行うための要素が2要素の数値という単純なデータであるがゆえではないかと予測します。例えば人間でも身長と体重データだけで性別や人種を分類できるかというと難しい気がしますが、他の要素も加えると精度は向上していきそうなことはイメージできます。
検証用データでも試してみる
最初にデータを準備するときに元のデータを、訓練用データと検証用データに分けておいたので、ここでは訓練用データによって最終的に得られた中心ベクトルを使用して検証用データのクラスタリングも行ってみます。
・実装例
plt.figure(figsize=(4,4)) R2=step1_kmeans(X_test[:,0],X_test[:,1],Mu) show_prm(X_test,R2,Mu,X_col)
・出力結果
ここでは中心ベクトル（Mu）は訓練用データにより得られた値を使用するので、Rの更新のみを一度行います。
また同様に正答率も出すと、「0.8267」となりました。
やはり正答率はそこまで高くない印象を受けますが、訓練用データの結果と概ね同程度であることから同種のデータに対しては同等の結果が得られるということが分かります。
だいぶサボってたので、pythonとブログの書き方を忘れた。
・今回のソース
python_dev/unsupervised_learning.ipynb at master · wantanblog/python_dev · GitHub
