実体と情報のはざま

何事にも囚われず。

散布図行列と主成分分析PCAをまとめ

 (*2017.06.14(18:30)ちょっとプログラムを修正しました。。)最近取り組んできた上記2つの手法は今後も頻繁に使うことになると思う。データサイエンティストを目指して修行を続けるなら。ただ、勉強を続ける中でデータサイエンティストとしての方向性が見えてきた。機械学習のような現実を近似するというか経験から物事を判別するようなアルゴリズムにはあまり興味がなくなってきた。データから実体の本質を見出す作業にこそやりがいを感じる。そういうわけで、今回のタイトルの手法はちゃんと押さえておきたいと思う。
 冗長なんだろうけど何をやっているか分かるプログラムに!
かなり長くなってしまった。。。データ行列のラベルに対するソートを四苦八苦しながらつくったもんで、素人っぽさがあふれだしてしまった。Pythonをはじめて2か月くらいたつだろうか。まだまだ、初心者の域は出ていないようだ。。
 内容はIrisデータセットの散布図行列と主成分分析。過去に両方ともやっているんだけれど、グラフに凡例を入れられなかったり、応用が利かない形になっていた部分をキレイにした。初心者にとっては簡単ではなかったけれど、冗長ながら納得のいくものができた。以前のものとぱっと見では同じかもしれないが、違いはディテールに表れている。カラーマップの使い方やなんかも作成過程でしっかり身についた。
まずは、散布図行列。
f:id:myuteru:20170614010721p:plain
お次が、主成分分析。
f:id:myuteru:20170614010740p:plain

 参考文献は相変わらずコレ。本当にお世話になっている。この本の著者くらいシンプルにコードが書けたらどんなに楽だろうと思う。。

 もう一つおすすめの本はコレ。基本的なことがしっかり盛りだくさん書かれていている。今まで紹介してこなかったけど、実はPythonをうちのパソコンにインストールする前から読み始めていた。はじめて買ったPython本。まじめな一冊。

 あと、プログラミングでくじけそうになった時の息抜きに読んですごく元気になれた本がコレ。すっごく前向きになれる。朝、通勤電車の中で読むと効果的!

 プログラムはコレ。長文注意!*2017.06.14(18:30)ちょっと修正しました。。
開発環境:Spyder(Python 3.6)

import numpy as np
import matplotlib.pyplot as plt

dat=np.genfromtxt('Iris.csv',delimiter=',',dtype=str)#まずは全部読み出し
X=np.copy(dat[1:,1:5])#数値データ行列を抜き出し
X=np.asfarray(X)#型変換str→float64
Y=np.copy(dat[1:,5])#ラベルの抜き出し
Z=np.copy(dat[0,1:5])#特徴量の名前

def sc_matrix(X,Y,Z):#散布図行列
    n=np.shape(X)
    L=np.zeros(n[0],dtype=int)#ラベルの数値化用
    c_name=np.empty(n[0],dtype=Y.dtype)#クラス名の取得用
    k,p=1,False
    for i in range(0,n[0]):
            if p==True:
                c_name[k]=Y[i-1]            
                k=k+1
            p=False
            for j in range(0,n[0]):
                if Y[i]==Y[j] and L[j]==0:
                    L[j]=k
                    p=True
    if p==True:
        c_name[k]=Y[i]
    Ln,k=L.max(),0#Ln:ラベルの種類数
    L_sort=np.zeros(n[0],dtype=int)#順番をソート
    cn=np.zeros(n[0],dtype=int)#クラス数のカウント
    for j in range(1,Ln+1):
        for i in range(0,n[0]):
            if L[i]==j:
                L_sort[k]=i
                k=k+1
        cn[j]=k
    X_sort=np.zeros((n[0],n[1]))
    Y_sort=np.copy(Y)
    for k in range(0,n[0]):
        X_sort[k,:]=np.copy(X[L_sort[k],:])
        Y_sort[k]=np.copy(Y[L_sort[k]])
    plt.figure(figsize=(14,12))#プロット開始
    k=0
    for j in range(0,n[1]):
       for i in range(0,n[1]):
           k=k+1
           plt.subplot(n[1],n[1],k)
           for m in range(0,Ln):
               plt.scatter(X_sort[cn[m]:cn[m+1],i],X_sort[cn[m]:cn[m+1],j],
                           s=15,alpha=0.8,c=plt.cm.cool((m+1)/Ln),
                           label=c_name[m+1])
           if j==n[1]-1:plt.xlabel(Z[i])#cmap=plt.cm.rainbow,c=m*col[cn[m]:cn[m+1]]
           if i==0 :plt.ylabel(Z[j])
           plt.legend(loc='best',fontsize='x-small',markerscale=0.5)
    plt.subplots_adjust(wspace=0.3,hspace=0.3)
    return

def pca(X,Y,Z):#主成分分析
    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]])#PC1
        G[i,1]=np.dot(X[i,:],ei_vec[:,a[1]])#PC2
    L=np.zeros(n[0],dtype=int)#ラベルの数値化用
    c_name=np.empty(n[0],dtype=Y.dtype)#クラス名の取得用
    k,p=1,False
    for i in range(0,n[0]):
            if p==True:
                c_name[k]=Y[i-1]            
                k=k+1
            p=False
            for j in range(0,n[0]):
                if Y[i]==Y[j] and L[j]==0:
                    L[j]=k
                    p=True
    if p==True:
        c_name[k]=Y[i]
    Ln,k=L.max(),0#ラベルの種類数
    L_sort=np.zeros(n[0],dtype=int)#順番をソート
    cn=np.zeros(n[0],dtype=int)#クラス数のカウント
    for j in range(1,Ln+1):
        for i in range(0,n[0]):
            if L[i]==j:
                L_sort[k]=i
                k=k+1
        cn[j]=k
    G_sort=np.zeros((n[0],2))
    Y_sort=np.copy(Y)
    for k in range(0,n[0]):
        G_sort[k,:]=np.copy(G[L_sort[k],:])
        Y_sort[k]=np.copy(Y[L_sort[k]])
    plt.figure(figsize=(15,12))#プロット開始
    pcax=np.empty(n[1],dtype=Y.dtype)
    q=np.arange(n[1])
    for i in range(0,n[1]):pcax[i]='PC'+str(i+1)
    plt.subplot(2,2,1)
    for m in range(0,Ln):
        plt.scatter(G_sort[cn[m]:cn[m+1],0],G_sort[cn[m]:cn[m+1],1],
                    s=30,alpha=0.8,c=plt.cm.cool((m+1)/Ln),
                    label=c_name[m+1])
    plt.legend(loc='best',fontsize='small',markerscale=0.5)    
    plt.xlabel("PC1",fontsize=12)
    plt.ylabel("PC2",fontsize=12)
    plt.subplot(2,2,2)
    plt.bar(q,ver,tick_label=pcax, align="center")
    plt.xlabel("principal components",fontsize=12)
    plt.ylabel("explained variance ratio",fontsize=12)
    plt.subplot(2,2,3)
    plt.bar(q,ei_vec[:,a[0]],tick_label=Z,align="center")
    plt.xlabel("PC1",fontsize=12)
    plt.ylabel("eigen vector",fontsize=12)
    plt.subplot(2,2,4)
    plt.bar(q,ei_vec[:,a[1]],tick_label=Z,align="center")
    plt.xlabel("PC2",fontsize=12)
    plt.ylabel("eigen vector",fontsize=12)
    plt.subplots_adjust(wspace=0.2,hspace=0.2)
    return

sc_matrix(X,Y,Z)
pca(X,Y,Z)

次は、少し毛色の違う手法である決定木とかの勉強をしたいな~!
今日はここまで!