LoginSignup
7
9

Delete article

Deleted articles cannot be recovered.

Draft of this article would be also deleted.

Are you sure you want to delete this article?

More than 5 years have passed since last update.

@MuAuan(ムー ウワン)

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

Last updated at Posted at 2019-02-26

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

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

やったこと

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

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

データは偏差値にしたものを利用しました。
普通にDistorsionは以下のようになっています。だからクラスタは３クラスか４クラス
keiba3__Distortion.jpg
そして、PCAの第一主成分と第二主成分面でのプロットは以下のとおり
※上記の回転を第一主成分vs第二主成分面でとめたものとみることもできる
keiba3_std2_PCA12_plotSn3.jpg

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

そして、第一主成分、第二主成分、そして第三主成分の寄与は以下の図のとおり
keiba3_contribution_PCA12_plotSn3.jpg
今回は、優勝、二着、三着のデータ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()
7
9
0

Register as a new user and use Qiita more conveniently

  1. You get articles that match your needs
  2. You can efficiently read back useful information
  3. You can use dark theme
What you can do with signing up
Sign upLogin
7
9

Delete article

Deleted articles cannot be recovered.

Draft of this article would be also deleted.

Are you sure you want to delete this article?