python を用いた主成分分析(PCA) のかんたんスクラッチ

#Ice Break

前回の記事では主成分分析(PCA)を題材に線形代数の総復習を行いました。

本稿ではpythonを用いたPCAの実装例について書きます。

sklearn使うとつまらないのでnumpy縛りで実装しました。

次回記事では今回作成したモジュールを用いて次元削減が回帰分析に対して与える影響について考察します。

#データセット準備
sklearnに付属しているBoston House Prices datasetを使用します。
まずhouse price 以外のデータを説明変数 X として使用するためにデータフレームからアレイとして取り出します。
また作法として X をカラムに沿って正規化しておきます。

from sklearn.datasets import load_boston
from sklearn.decomposition import PCA
import pandas as pd
import numpy as np

def create_X():
    boston = load_boston()
    data = pd.DataFrame(boston.data, columns = boston.feature_names)
    data = data.values
    X = np.delete(data,12,1)
    mean = np.mean(X,axis=0)
    std  = np.std(X,axis=0) 
    X = (X-mean)/std
    return X

#PCA
前回の記事で述べた通り、PCAでは$X^{T}X$に対する固有ベクトルを探索します。

そのためまず始めに$X^{T}X$を計算します。(1)

つぎにnumpy.linali.eigを用いて固有ベクトルを計算します。(2)

$D$は固有ベクトルを対応する固有値が大きい順に並べたものです。任意の次元に圧縮するためn_componensを変数としています。また以下で数式と同様に計算を進めるために一度転置を挟んでいますが、実装上は必要ないです。(3)

最後に$D^{T}X$を計算することで次元を削減することが出来ます。(4)

def PCA_Sc(X,n_components):
    columns = X.shape[1]
    #　1 calucurate X_TX
    X_TrX = np.dot(X.T,X)
    # 2 clucurate eigenvalue & Eigenvector
    w,v = np.linalg.eig(X_TrX)
    # 3 create D
    D_list = [v[:,i] for i in range(n_components)]
    D_array = np.array(D_list)
    D_reshape = np.reshape(D_array,(n_components,columns))
    D = D_reshape.T
    # 4 colucurate D_TX
    X_transform = np.dot(D.T,X.T)
    return X_transform

#sklearn との比較
sklearnで生成した主成分$D$との比較を行いました。結果は以下になります。

この図から分かることはsklearnでは第１主成分の最初の値が正になるように符号を反転させる機能が内蔵されている事です。（このことは他の方も報告していたと思います。）

スクラッチで書くとこのような発見があって楽しいですね。
定義参照したほうが早いですが…