【Python】 主成分分析（PCA）
導入：データの次元削減
主成分分析の目的
主成分分析（principal component analysis，PCA）とは一体何でしょうか？PCA には大きく分けて二つの目的があります。（表していることは同じですが）
- データの特徴をよく表す新しい指標を見つける
- 大きな次元のデータを小さな次元に落とす
データの次元が大きすぎて、計算にやたら時間がかかる場合とかによく使ったりします。
多次元のデータを解析する
世の中に存在するデータは多次元のものがほとんどです。各次元の値は特徴量と呼ばれたりします（呼ばせてください）。例えば、「人間」というデータは、“身長”、“体重”、“年齢”、“名前”などの特徴量を持っています。これが学生であれば、各科目のテストの点数なども含まれるかもしれません。
もし、この特徴量の中に“目の数”が含まれていたらどうでしょうか。これはほとんど意味がないですよね。このように、特徴量のなかにはそのデータにとって重要なものや、あまり重要でないものが含まれます。しかも、どの特徴量がどれくらい重要かなどを人力で決定するのは不可能ですよね。
データをよく表す指標（主成分）を見つける
そこで、統計的にデータを解析して、データをよく表す指標（主成分）を見つけるアルゴリズムが主成分分析です。「統計的にデータを解析する」とは、たくさんデータを集めて、それらを数学的に解析するという意味です（適当）。
より具体的なイメージをつかむために、データが二次元の場合の例をみてみましょう。次の図では、青い各点が一つのデータを表しています。
この図を見ると、赤い矢印方向の軸がこのデータをよく表しているのが分かりますね。この赤い矢印を見つけるのが主成分分析です。赤色の軸にデータを写像することで、次元を２から１に減らすことができますね。
主成分分析のやり方
上の例は次元を２から１に落とすものですが、主成分分析は写像したとき（次元を落としたとき）になるべくデータの分散が大きくなるような軸を探すタスクです。主成分の軸と元データの二乗誤差が最小になるようにすると説明されていることもありますが、これはデータの平均が原点と同じ場合であり、このとき分散を最大化することと等価になります。なので、主成分分析は- データを、平均が原点になるようにシフトする
- そのデータの二乗誤差が最小となるような軸を探す
という流れで行うことが多いと思います。（平均を原点にしておかないと全然意味の無い軸が出てきてしまうことがある。）
では具体的に、主成分の軸はどうやって計算するのかという話ですが、実は簡単で、各データが列に格納されているような行列 M (∈ dim * N)を考えたとき、M の共分散行列の固有値と固有ベクトルを計算し、固有値の大きい順に第一主成分、第二主成分…となっており、対応する固有ベクトルが主成分軸です。 なんでそれで主成分が計算できるのか気になりますよね。式とか書くのが面倒くさいので、適当に調べてください。
適当に調べて出てきた資料ですが、割と分かりやすいと思います。
http://pub.maruzen.co.jp/realize/search/sample/RLB220479.pdf
次元が大きすぎる場合
元データの次元がデータ数に比べてとても大きい場合があります。例えば画像データなどはピクセル数が次元数となります。こういった場合、共分散行列が次元数×次元数の行列となり、メモリも足りませんし、計算時間も大変なことになってしまいます。そこで転置したデータ行列の共分散行列を計算する方法があります。これだと共分散行列はデータ数×データ数となり、なんとかなりそうですね。
下の実装では、データ数と次元数によってこれらを使い分けているので、参照してください。
実装
主成分分析の実装です。データを何次元に落としたいか（第何主成分まで欲しいか）を引数として与えることができます。返り値は、主成分軸の基底です。# -*- coding: utf-8 -*- from matplotlib import pylab as pl from numpy.random import randn from scipy import array, ceil, dot, float64, floor, matrix, sqrt, zeros from scipy.linalg import cholesky, eig, norm, solve # ================= # 主成分分析 (PCA) # ================= """ Maximum Variance criterion """ def pca(data, base_num = 1): N, dim = data.shape data_m = data.mean(0) data_new = data - data_m ### データ数 > 次元数 if N > dim: ### データ行列の共分散行列 cov_mat = dot(data_new.T, data_new) / float(N) ### 固有値・固有ベクトルを計算 l, vm = eig(cov_mat) ### 固有値が大きい順に並び替え axis = vm[:, l.argsort()[- min(base_num, dim) :][:: -1]].T ### 次元数 > データ数 else: base_num = min(base_num, N) cov_mat = dot(data_new, data_new.T) / float(N) l, v = eig(cov_mat) ### 固有値と固有ベクトルを並び替え idx = l.argsort()[::-1] l = l[idx] v = vm[:, idx] ### 固有ベクトルを変換 vm = dot(data_m.T, v[:, : base_num]) ### （主成分の）基底を計算 axis = zeros([base_num, dim], dtype = float64) for ii in range(base_num): if l[ii] <= 0: break axis[ii] = vm[:, ii] / norm(vm[:, ii]) return axis # ======== # テスト # ======== from numpy.random import multivariate_normal def test(): data = multivariate_normal([0, 0], [[1, 2], [2, 5]], 100) ### PCA pc_base = pca(data, base_num = 1)[0] ### Plotting fig = pl.figure() fig.add_subplot(1,1,1) pl.axvline(x=0, color = "#000000") pl.axhline(y=0, color = "#000000") ### Plot data pl.scatter(data[:, 0], data[:, 1]) ### Draw the 1st principal axis pc_line = array([-3., 3.]) * (pc_base[1] / pc_base[0]) pl.arrow(0, 0, -pc_base[0] * 2, -pc_base[1] * 2, fc = "r", width = 0.15, head_width = 0.45) pl.plot([-3, 3], pc_line, "r") ### Settings pl.xticks(size = 15) pl.yticks(size = 15) pl.xlim([-3, 3]) pl.tight_layout() pl.show() if __name__ == "__main__": test()
テストコードが無駄に長くなっていますが、PCA の関数は非常に簡単に書けていることが分かります。
