【機械学習入門】「PCA+K-meansクラスタリング」を３D見える化♬

頭の中ではきっとこんなイメージ；PCA＋K-meansクラスタリング

Qiitaへのアップの都合上、若干サイズや分解能を犠牲にしたが、うまく表示できました。
つまり、第一主成分が一番広がりが大きく（分散が大きく）その広がりの方向に延びており（固有値問題の解）、二番目が第二主成分（固有値問題の二番目の解）。。。となっており、今回のように各要素を偏差値で置き換えてPCA分析するのが正解なのが分かる。
【参考】
・[Python][Scikit-learn]主成分分析を用いた次元削減、主成分ベクトルを用いた予測と線形回帰による予測の比較

やったこと

・競馬騎手データをとりあえず３次元データのみを利用してクラスタリングしてみる
・３ｄプロットしてみる

・競馬騎手データを３次元のみのデータを利用してクラスタリングしてみる

データは偏差値にしたものを利用しました。
普通にDistorsionは以下のようになっています。だからクラスタは３クラスか４クラス

そして、PCAの第一主成分と第二主成分面でのプロットは以下のとおり
※上記の回転を第一主成分vs第二主成分面でとめたものとみることもできる

相関は以下のとおり
※これも回転して1stvs2nd、2ndvs3rd, そして3rdvs1st面で止めた絵になっているはずです

そして、第一主成分、第二主成分、そして第三主成分の寄与は以下の図のとおり

今回は、優勝、二着、三着のデータ120個を使っているので、上記のように相関が高く、第一主成分の寄与が８８％以上と高く、第二と第三が６％、５％と小さいです。
ちなみに、第一、二、三主成分の軸方向のベクトル座標は以下のようになっています。
※固有値だからお互い直交（ベクトルの内積０）してます

v_pca0=[ 0.57958236  0.5776769   0.57478143]
v_pca1=[ 0.25818473  0.53882344 -0.80187901]
v_pca2=[ 0.77293269 -0.61315472 -0.1631452 ]

・３ｄプロットしてみる

今回は、少しずるしてもともとの座標で描画しています。
ということで、肝心な部分のコードは以下のとおり
※コード全体はおまけに記載します
まず、入力データに対して主成分分析をpcaで実施します。
それぞれの第一主成分などはpca.components_[0]などで求められます。

#主成分分析の実行
pca = PCA()
pca.fit(df.iloc[:, 1:])
PCA(copy=True, n_components=None, whiten=False)
v_pca0 = pca.components_[0]
v_pca1 = pca.components_[1]
v_pca2 = pca.components_[2]

出力してみると、元の座標軸での座標で示されています。
そこから、X座標などを分割します。
※ベクトルは本来原点からの座標になっているので、今回の偏差値に合わせて原点を（50，50，50）にずらします
第三主成分までのベクトルを引きたいのでそれぞれの座標（X,Y,Z）として以下のように求めます

print("v_pca0={},v_pca1={},v_pca2={}".format(v_pca0,v_pca1,v_pca2))
X_pca1=[50,50+30*v_pca0[0]]
Y_pca1=[50,50+30*v_pca0[1]]
Z_pca1=[50,50+30*v_pca0[2]]
X_pca2=[50,50+30*v_pca1[0]]
Y_pca2=[50,50+30*v_pca1[1]]
Z_pca2=[50,50+30*v_pca1[2]]
X_pca3=[50,50+30*v_pca2[0]]
Y_pca3=[50,50+30*v_pca2[1]]
Z_pca3=[50,50+30*v_pca2[2]]

ラベルの表示もしたいので、以下のように定義します。

# 分類結果のラベルを取得する
labels = np.unique(kmeans_model.labels_)

そして、これらを使ってangle毎にプロットを表示します。

def show_with_angle(angle):
    fig = plt.figure()
    ax = fig.add_subplot(111, projection='3d')

カメラ視点を以下のとおり定義します。
【参考】
・how to set “camera position” for 3d plots using python/matplotlib?

    ax.view_init(elev = 10., azim = angle) #視点

以下で、求まったクラスタの中心近くにlabelsで定義した名前（０，１，２）をtext記載します。

    for kx, ky,kz, name in zip(kmeans_model.cluster_centers_[:,0], kmeans_model.cluster_centers_[:,1], kmeans_model.cluster_centers_[:,2], labels):
        ax.text(kx, ky,kz, name, alpha=0.8, size=20)

そして上記で定義した第一成分などの軸を記載します。

    ax.plot(X_pca1,Y_pca1,Z_pca1, c='b', marker='o', alpha=0.5, label='1st pricipal')  
    ax.plot(X_pca2,Y_pca2,Z_pca2, c='r', marker='o', alpha=0.5, label='2nd principal') 
    ax.plot(X_pca3,Y_pca3,Z_pca3, c='g', marker='o', alpha=0.5, label='3rd principal') 

以下で本来のデータを分類されたcolorsの色でプロットします。

    ax.scatter(df.iloc[:, 1], df.iloc[:, 2], df.iloc[:, 3],  c=colors, marker='o', alpha=1)

次にクラスタの中心に大き目（s=100）なプロットをします。

    ax.scatter(kmeans_model.cluster_centers_[:,0], kmeans_model.cluster_centers_[:,1], kmeans_model.cluster_centers_[:,2], c = "b", marker = "*", s = 100)

最後に軸ラベルと凡例（凡例は上記のプロットのlabelで定義したもの）を記載します。

    ax.set_xlabel('1st')
    ax.set_ylabel('2nd')
    ax.set_zlabel('3rd')
    ax.legend()
    plt.pause(0.01)
    plt.savefig('k-means/keiba3/pca3d/keiba3_PCA3d_angle_'+str(angle)+'.jpg')
    plt.close()

角度は360度を1度ずつ計算しました。
容量の関係で上記のGifアニメーションは5度ずつのデータを表示しています。

for angle in range(0, 360, 1):
    show_with_angle(angle)

まとめ

・PCA＋K-meansクラスタリングを３D見える化してみた
・思い通りのGifアニメーションとなった

・さらなら高次元の見える化に挑戦したい

おまけ

import pandas as pd # データフレームワーク処理のライブラリをインポート
import matplotlib.pyplot as plt
from pandas.tools import plotting # 高度なプロットを行うツールのインポート
from sklearn.cluster import KMeans # K-means クラスタリングをおこなう
from sklearn.decomposition import PCA #主成分分析器
from mpl_toolkits.mplot3d import Axes3D
import numpy as np


df = pd.read_csv("keiba100_std3.csv", sep=',', na_values=".") # データの読み込み
df.iloc[:, 1:].head() #解析に使うデータは２列目以降
print(df.iloc[:, 1:])
X=df.iloc[:, 1:]

plotting.scatter_matrix(df[df.columns[1:]], figsize=(6,6), alpha=0.8, diagonal='kde')   #全体像を眺める
plt.savefig('k-means/keiba3/keiba3_scatter_plot4.jpg')
plt.pause(1)
plt.close()

distortions = []
distortions1 = []
for i  in range(1,21):                # 1~20クラスタまで一気に計算 
    km = KMeans(n_clusters=i,
                init='k-means++',     # k-means++法によりクラスタ中心を選択
                n_init=10,
                max_iter=300,
                random_state=0)
    km.fit(df.iloc[:, 1:])                # クラスタリングの計算を実行
    distortions.append(km.inertia_)       # km.fitするとkm.inertia_が得られる
    UF = km.inertia_ + i*np.log(20)*5e2   # km.inertia_ + kln(size)
    distortions1.append(UF)   

fig=plt.figure(figsize=(12, 10))
ax1 = fig.add_subplot(311)
ax2 = fig.add_subplot(312)
ax3 = fig.add_subplot(313)
ax1.plot(range(1,21),distortions,marker='o')
ax1.set_xlabel('Number of clusters')
ax1.set_ylabel('Distortion')
ax2.plot(range(1,21),distortions,marker='o')
ax2.set_xlabel('Number of clusters')
ax2.set_ylabel('Distortion')
ax2.set_yscale('log')
ax3.plot(range(1,21),distortions1,marker='o')
ax3.set_xlabel('Number of clusters')
ax3.set_ylabel('Distortion+klog')
ax3.set_yscale('log')
plt.pause(1)
plt.savefig('k-means/keiba3/keiba3__Distortion.jpg')
plt.close()

s=3
# この例では s個のグループに分割 (メルセンヌツイスターの乱数の種を 10 とする)
kmeans_model = KMeans(n_clusters=s, random_state=10).fit(df.iloc[:, 1:])

# 分類結果のラベルを取得する
labels = kmeans_model.labels_
print(labels)     ##################  print(labels) ###################
# それぞれに与える色を決める。
color_codes = {0:'#00FF00', 1:'#FF0000', 2:'#0000FF', 3:'#00FFFF' , 4:'#007777', 5:'#f0FF00', 6:'#FF0600', 7:'#0070FF', 8:'#08FFFF' , 9:'#077777',10:'#177777', 11:'#277777', 12:'#377777', 13:'#477777' , 14:'#577777'}
# サンプル毎に色を与える。
colors = [color_codes[x] for x in labels]

fig, axes = plt.subplots(nrows=3, ncols=3, figsize=(18, 18),sharex=False,sharey=False)
for i in range(1,4):
    for j in range(1,4):
        ax = axes[i-1,j-1]
        x_data=df[df.columns[i]]
        y_data=df[df.columns[j]]
        ax.scatter(x_data,y_data,c=colors)
        cor=np.corrcoef(x_data,y_data)[0, 1]
        ax.set_title("{:.3f}".format(cor),fontsize=30)
plt.savefig('k-means/keiba3/keiba3_correlation.jpg')
plt.close()

#主成分分析の実行
pca = PCA()
pca.fit(df.iloc[:, 1:])
PCA(copy=True, n_components=None, whiten=False)
v_pca0 = pca.components_[0]
v_pca1 = pca.components_[1]
v_pca2 = pca.components_[2]

print("v_pca0={},v_pca1={},v_pca2={}".format(v_pca0,v_pca1,v_pca2))
X_pca1=[50,50+30*v_pca0[0]]
Y_pca1=[50,50+30*v_pca0[1]]
Z_pca1=[50,50+30*v_pca0[2]]
X_pca2=[50,50+30*v_pca1[0]]
Y_pca2=[50,50+30*v_pca1[1]]
Z_pca2=[50,50+30*v_pca1[2]]
X_pca3=[50,50+30*v_pca2[0]]
Y_pca3=[50,50+30*v_pca2[1]]
Z_pca3=[50,50+30*v_pca2[2]]

# 分類結果のラベルを取得する
labels = np.unique(kmeans_model.labels_)

feature = pca.transform(df.iloc[:, 1:])
x,y,z = feature[:, 0], feature[:, 1], feature[:, 2]
X=np.c_[x,y,z]

def show_with_angle(angle):
    fig = plt.figure()
    ax = fig.add_subplot(111, projection='3d')
    ax.view_init(elev = 10., azim = angle)
    for kx, ky,kz, name in zip(kmeans_model.cluster_centers_[:,0], kmeans_model.cluster_centers_[:,1], kmeans_model.cluster_centers_[:,2], labels):
        ax.text(kx, ky,kz, name, alpha=0.8, size=20)
    ax.plot(X_pca1,Y_pca1,Z_pca1, c='b', marker='o', alpha=0.5, label='1st pricipal')  
    ax.plot(X_pca2,Y_pca2,Z_pca2, c='r', marker='o', alpha=0.5, label='2nd principal') 
    ax.plot(X_pca3,Y_pca3,Z_pca3, c='g', marker='o', alpha=0.5, label='3rd principal') 
    ax.scatter(df.iloc[:, 1], df.iloc[:, 2], df.iloc[:, 3],  c=colors, marker='o', alpha=1)
    ax.scatter(kmeans_model.cluster_centers_[:,0], kmeans_model.cluster_centers_[:,1], kmeans_model.cluster_centers_[:,2], c = "b", marker = "*", s = 100)
    ax.set_xlabel('1st')
    ax.set_ylabel('2nd')
    ax.set_zlabel('3rd')
    ax.legend()
    plt.pause(0.01)
    plt.savefig('k-means/keiba3/pca3d/keiba3_PCA3d_angle_'+str(angle)+'.jpg')
    plt.close()

for angle in range(0, 360, 1):
    show_with_angle(angle)

print("pca.explained_variance_ratio_={}".format(pca.explained_variance_ratio_))
plt.bar([1, 2,3], pca.explained_variance_ratio_, align = "center")
plt.title("pca.explained_variance_ratio_={}".format(pca.explained_variance_ratio_))
plt.xlabel("components")
plt.ylabel("contribution")
plt.savefig('k-means/keiba3/keiba3_contribution_PCA12_plotSn'+str(s)+'.jpg')
plt.pause(1)
plt.close()

# データを主成分空間に写像 = 次元圧縮
feature = pca.transform(df.iloc[:, 1:])
x,y = feature[:, 0], feature[:, 1]
X=np.c_[x,y]
kmeans_model = KMeans(n_clusters=s, random_state=10).fit(X) #PCA平面（ｘ、ｙ）で再度クラスタリング

plt.figure(figsize=(6, 6))
for kx, ky, name in zip(x, y, df.iloc[:, 0]):
    plt.text(kx, ky, name, alpha=0.8, size=10)
plt.scatter(x, y, alpha=0.8, color=colors)
plt.scatter(kmeans_model.cluster_centers_[:,0], kmeans_model.cluster_centers_[:,1], c ="b" , marker = "*", s = 200)
plt.title("Principal Component Analysis")
plt.xlabel("The first principal component score")
plt.ylabel("The second principal component score")
plt.savefig('k-means/keiba3/keiba3_std2_PCA12_plotSn'+str(s)+'.jpg')
plt.pause(1)
plt.close()