実体と情報のはざま

何事にも囚われず。

PCA(主成分分析)を手中に収める!

 機械学習を学んでいて、やたらと出てくるのがこのPCAというやつ。どの文献でもさらっと事前準備みたいな扱い。「まずはPCAで次元を削減してから・・・」なんてよく出てくる。ありきたりな手法のようですな。しかし、データサイエンス初心者の私としては、なんであれ一通りは理解しながら前へ進みたい!今回は、このPCAの学習と実装について書き残すことにする。
ライブラリを使えば簡単ですが!
 Pythonではおなじみのsckit-learnには簡単に実装可能なライブラリが用意されている。しかし、いつものことながら理論的な部分もしっかり学びたいので、ライブラリは使わないでやりたい。ただ、ここから少し深みにはまった。PCAの理論を解説している文献はいくつもネットに転がっているものの、理論の理解から実装まで導いてくれるものはなかった。理論の部分はネットの参考文献でよさそうなものを数件選んで勉強した。実装の部分は、自分の本棚からPython機械学習プログラミング 達人データサイエンティストによる理論と実践 (impress top gear)を引っ張り出してきて参考にした。参考にはしたものの出来上がったものはちょっと違うものになっている。というのも、自分が理解した理論の流れに沿ってプログラミングをしていくと、どうしても参考文献とは異なってしまうから。ここに自己流が入ることが不安でもあり楽しいところでもある。
できた!けどあってるのか!?
 データセットはデータサイエンスでは有名なIrisとwineを使った。前者は以前kaggleから入手したが、wineはUCI Machine Learning Repositoryから入手した。このサイトも有名。実はIrisデータセットはここにもある。できたPCAのプログラムはこんな感じ。
開発環境:Spyder(Python 3.6)

import numpy as np
import matplotlib.pyplot as plt

Iris=np.genfromtxt('Iris.csv',delimiter=',',dtype=str)#まずは全部読み出し
X=np.copy(Iris[1:,1:5])#データ行列を抜き出し
X=np.asfarray(X)#型変換str→float64

def PCA(X):
    n=X.shape#データ行列数(行,列)を取得
    for j in range(0,n[1]):
        X[:,j]=(X[:,j]-X[:,j].mean())/X[:,j].std()#データの標準化
    cov_X=np.cov(X.T)#共分散行列
    ei_val,ei_vec=np.linalg.eigh(cov_X)#固有値と固有ベクトル
    tot=np.sum(ei_val)#固有値の和
    s_ei_val=sorted(ei_val,reverse=True)#固有値の並べ替え
    ver=s_ei_val/tot#分散説明率
    a=np.zeros(n[1],dtype=np.int32)
    b=np.copy(np.abs(ei_val))
    for i in range(n[1]):#固有値の大きい順に並べ替える
        a[i]=np.argmax(b)#bが最大となるインデックスを格納
        b[a[i]]=0#値を取り出したらゼロにする
    g=np.zeros((n[0],2))
    for i in range(0,n[0]):
        g[i,0]=np.dot(X[i,:],ei_vec[:,a[0]])
        g[i,1]=np.dot(X[i,:],ei_vec[:,a[1]])
    return ver,ei_vec,a,g

ver,ei_vec,a,g=PCA(X)
n=X.shape
##############以下、グラフ化作業
#rmd=np.dot(np.dot(ei_vec[:,a[1]],cov_X),ei_vec[:,a[1]])
plt.figure(0)
plt.bar(np.arange(n[1]),ver[0:n[1]])
plt.xlabel("principal components")
plt.ylabel("explained variance ratio")
plt.figure(1)
plt.bar(np.arange(n[1]),ei_vec[:,a[0]])
plt.xlabel("PC1")
plt.ylabel("Y")
plt.figure(2)
plt.bar(np.arange(n[1]),ei_vec[:,a[1]])
plt.xlabel("PC2")
plt.ylabel("Y")
plt.figure(3)
for i in range(0,n[0]):
    if Iris[i+1,5]=='Iris-setosa':
        plt.scatter(g[i,0],g[i,1],color='b',alpha=0.7)
    elif Iris[i+1,5]=='Iris-versicolor':
        plt.scatter(g[i,0],g[i,1],color='r',alpha=0.7)
    elif Iris[i+1,5]=='Iris-virginica':
        plt.scatter(g[i,0],g[i,1],color='g',alpha=0.7)
plt.xlabel("PC1")
plt.ylabel("PC2")
plt.legend()

 結果はこんな感じ。次元削減というよりは教師なし分類として見ると、グラフが単色で見分けがつくのは青いやつだけかな。ところで、他のサイトで同じことをやっているはずのプロットと若干異なる。なんで?よく見ると、ネットに転がっている同様の結果どうしも若干違っている!そんなものなのか?今はいいや。いつか振り返ってみよう。プログラム中には分散説明率のグラフも入っているので忘備録として解説を入れたいが説明が長くなりそうなので省略。
f:id:myuteru:20170607232450p:plain
 つぎに、wineデータセットを使った場合の結果。これも教師なし分類としてみると単色で見分けがつくのは緑だけかな。ちなみに、上のプログラムのdef以外の部分をwineデータセットに合わせてちょいといじるだけでできる。
f:id:myuteru:20170607233029p:plain
 今回は使い古されてなお現役で機械学習の前座で重宝されているPCAを実装した。この手法自体をちょっと軽んじてたけど、実際に四苦八苦しながらやってみると、初めに考えた人のすごさや奥深さを知ることになり大変に勉強になった。。
 おっと、日付が変わってしまう!今日はここまで!