【python】sklearnのPCAでloading（主成分負荷量）を計算する

　PCA（主成分分析）のloading*1がほしいときがあります。

　sklearnでは一発では出ません。

　ドキュメントはここ。
sklearn.decomposition.PCA — scikit-learn 0.21.2 documentation

PCA.components_は確かにあるけど・・・
loadingを計算しよう
罠だった
共分散行列のときはどうするのか
loadingを使うと何が良いのか

PCA.components_は確かにあるけど・・・

　結論から先に言うと、PCA.components_はノルム1の固有ベクトルです。ノルムを測ってみましょう。

>>> import numpy as np
>>> from sklearn.datasets import load_iris
>>> from sklearn.decomposition import PCA
>>> iris = load_iris()
>>> pca = PCA(n_components=2)
>>> pca.fit(iris.data)
PCA(copy=True, iterated_power='auto', n_components=2, random_state=None,
  svd_solver='auto', tol=0.0, whiten=False)
>>> pca.components_
array([[ 0.36158968, -0.08226889,  0.85657211,  0.35884393],
       [ 0.65653988,  0.72971237, -0.1757674 , -0.07470647]])
>>> np.linalg.norm(pca.components_, axis=1)
array([1., 1.])

　まあ、loadingも固有ベクトルには違いないのですが、ノルムを整えてやる必要があります。

loadingを計算しよう

　教科書などによく書いてあることですが、第 $i$ 主成分に対応する元の変数のloadingは次の式で出せます。
$loading = \sqrt(\lambda_i) eigenvector$
　 $\lambda_i$ は固有値。 eigenvectorは固有ベクトルで、元の変数の数だけ次元がありますから、これで良いわけです（雑な説明ですが・・・）。

　ということで、pythonで同様にやってみましょう。固有値はexplained_varianceに入っています。

>>> pca.components_*np.c_[np.sqrt(pca.explained_variance_)] # 縦ベクトルに変換する必要がある
array([[ 0.74322652, -0.16909891,  1.76063406,  0.73758279],
       [ 0.32313741,  0.35915163, -0.08650963, -0.03676921]])

　できました。これがloadingです。・・・と思ったけど、1を超えちゃってますね。なんてこった。

罠だった

　固有値は分散なので、データのスケールに依存します。

　とりあえず入力をスケーリングしてみよう。上の式は相関行列から行くときのものでした。なのでこれで平気なはず。

>>> from sklearn.preprocessing import StandardScaler as SS
>>> ss = SS()
>>> data = ss.fit_transform(iris.data)
>>> pca.fit(data)
PCA(copy=True, iterated_power='auto', n_components=2, random_state=None,
  svd_solver='auto', tol=0.0, whiten=False)
>>> pca.components_*np.c_[np.sqrt(pca.explained_variance_)]
array([[ 0.89421016, -0.45081822,  0.99500666,  0.96822861],
       [ 0.35854928,  0.89132754,  0.02031465,  0.06299656]])

　1を超えなくなくてめでたし、ということよりも、数字が変わったことの方が問題で、これで本当に正しいのかという疑念が生じてきました。

　確認のために元の特徴と主成分の相関係数を直接測ってみます。

>>> X = pca.fit_transform(data)
>>> np.corrcoef(np.hstack([iris.data, X]), rowvar=False)
array([[ 1.00000000e+00, -1.09369250e-01,  8.71754157e-01,  8.17953633e-01,  8.91224479e-01,  3.57352114e-01],
       [-1.09369250e-01,  1.00000000e+00, -4.20516096e-01,  -3.56544090e-01, -4.49312976e-01,  8.88351481e-01],
       [ 8.71754157e-01, -4.20516096e-01,  1.00000000e+00,  9.62757097e-01,  9.91684422e-01,  2.02468206e-02],
       [ 8.17953633e-01, -3.56544090e-01,  9.62757097e-01,  1.00000000e+00,  9.64995787e-01,  6.27862218e-02],
       [ 8.91224479e-01, -4.49312976e-01,  9.91684422e-01,  9.64995787e-01,  1.00000000e+00,  2.08904471e-17],
       [ 3.57352114e-01,  8.88351481e-01,  2.02468206e-02,  6.27862218e-02,  2.08904471e-17,  1.00000000e+00]])

　下の二行の4列目までを見てください。微妙に誤差があるような気はしますが（小数点以下3桁でずれてきてるので微妙ってほど微妙でもないが）、たぶん同じ数字になっています。

　微妙な誤差については、丸め誤差などが累積した、実は計算間違ってる、といった理由が考えられます。前者ならまだ許せるけど、後者はやだな・・・。

共分散行列のときはどうするのか

　どうするんだろうね・・・。

　2019/06/22追記
　手順は増えますが、スケールを考慮すれば同様に行えるようです。

出典：http://manabukano.brilliant-future.net/document/text-PCA.pdf　p.10

loadingを使うと何が良いのか

　相関係数なので、「どれくらい効いてるか」がよくわかります。よくある0.3以下なら～とか0.7以上なら～という論法ができます。それだけといえばそれだけ。

　このように取扱が大変なので、固有ベクトルのまま解釈した方が楽かもという気がしてきました。各主成分の寄与率はexplained_variance_ratio_で得られる訳だし、寄与率の大きい軸の固有ベクトルの大きい次元を見ればどんな感じかはわかるし・・・。

　でもまあ、一応（入力をスケーリングすれば）大体出るということはわかったので、これでよしとします。

　共分散行列でやるときのやり方は、どなたか詳しい方に教えて頂けると助かります。

*1:主成分負荷量、あるいは因子負荷量とも（なぜか）言われますが、この記事ではloadingで通します。けっきょくヘタに和訳しないのがいちばんわかりやすい

静かなる名辞

pythonとプログラミングのこと