機械学習
線形代数
主成分分析

線形な手法からカーネル法へ(2. 主成分分析)

⚠️記事のアップデートの予定あり


はじめに

前回は線形回帰とカーネル回帰について説明しました。今回は主成分分析(PCA: Principal Component Analysis)について解説します。

今回は主成分分析について解説します。主成分分析は100年くらい前からある古典的な手法ですが、たとえば最近流行りのニューラルネットワークの中でも白色化という形で使われたりしています。主成分分析についてはいろいろ説明することがあるので、今回はなるべくミニマルな説明に限定し、残りの話は番外編として書きます。

それでは始めましょう。


注意

主成分分析は既に詳しい解説も多いです。Wikipediaは一度目を通す価値があるでしょう。

他にもよい記事[1]がたくさんあるので、自分で理解しやすい記事を参考にするとよいと思います。


主成分分析

主成分分析とは、変換後の各確率変数の分散が最大になるような直交変換を求めることです。主成分分析を行う場合には次の仮定が満たされているか確認する必要があります。


  1. 観測したい情報の持つ分散のほうがノイズに起因する分散よりも非常に大きい。

  2. 観測したい真の情報の次元は実際に観測された情報の次元よりも小さい。

主成分分析がうまくいくかどうかは仮定1が満たされているかどうかに依存しています。主成分分析は確率変数の分散が最大となるような方向を取り出す操作なので、ノイズが大きいとまったくデタラメな方向を抽出してしまい役に立ちません。また主成分分析はよく次元圧縮の手法のひとつとして紹介されますが、主成分分析の本質は「確率変数の無相関化」です。次元圧縮が可能かどうかは仮定2が満たされているかどうかに依存します。


注意:「確率変数が無相関であること」と「確率変数が独立であること」は必ずしも一致しません(参考)。


データセット

さて、この時点でもデータ分析の言葉に慣れていない人にとっては難しい言葉がいくつかあるように見えるかもしれませんが、具体例を挙げながら解説していきますのでご安心ください。

たとえば学校で生徒の学業の優秀さを測るために数学と英語のテストを実施したとします(本当はもっとたくさんの科目を実施するでしょうが、可視化しやすいように2次元に絞ります)。学校は生徒を成績によって序列づけることをしますので、このテストの結果に基づいて生徒に順位を割り当てようとしています。優劣の比較をするためには1次元に落とし込む必要があるので、これは「生徒の学業の優秀さは本当は1次元で表現できるはずだ」という仮定に基づく「2次元から1次元への次元圧縮」の問題です。次元圧縮(あるいは次元削減)とは、データの持つ本質的な情報を保ったまま多次元のデータからより小さい次元のデータへ変換することを言います。

test_scores.png

上のデータは「生徒の学業の優秀さを表す潜在的な方向ベクトル${\bf a}\,(\|{\bf a}\|=1,\,$図中の青いベクトル$)$が存在する。その係数はガウス分布$G(\mu, \sigma)$に従い、さらにテストの点数には観測誤差として2次元のガウスノイズが加わる」というモデルにしたがって生成したものです。つまり数学と英語の結果に意図的に正の相関を持たせた「数学の出来がいい子は英語の出来もいい」という倫理的な物議を醸しかねないデータを説明のために仕方なく生成したわけです。上の図には生徒たちの成績および平均点${\bf m}$と潜在的な方向ベクトル${\bf a}$もプロットしてあります(${\bf a}$は平均点の位置を始点とし、見やすいように拡大してあります)。この${\bf a}$が不明なときに${\bf a}$と${\bf m}$それぞれを復元することができれば主成分分析は成功というわけです。

学校で通常用いられる単純な合計点と主成分分析との比較がしやすいように、このデータは実際にありえるかありえないかギリギリのラインまで意図的に${\bf a}$を偏らせてあります。つまり数学の教師はけっこう厳しい先生なので下は0点付近まで取りえるようなテストになっていますが、英語の教師は慈悲深い先生なので50点以下はほぼいないだろうというテストになっています。生徒のテストの結果は次のコードで生成しました。


input

import numpy as np

from matplotlib import pyplot as plt
import seaborn as sns

# m個のn次元のデータを生成する(mxn行列)
# Xの行数はp(p < n)とし、Xは行ごとに主成分ベクトルを持つ。
# mu と sigma はp個ずつ指定し、0...p-1番目は主成分の平均と分散である。
# noise_mu と noise_sigma はn個ずつ指定し、0...n-1番目は観測軸に沿ったノイズの平均と分散である。
# こうして生成されたp本のペクトルをXで重みづけて線型結合し、
# ノイズを加えたものを実際に観測されたデータとする。
def make_dataset(X, mu, sigma, noise_mu, noise_sigma, m, n):
pc_coefs = np.zeros((m, X.shape[0]))
for i in range(X.shape[0]):
pc_coefs[:,i] = np.random.normal(mu[i], sigma[i], m)
score = pc_coefs.dot(X)

noise = np.zeros((m, n))
for i in range(n):
noise[:,i] = np.random.normal(noise_mu[i], noise_sigma[i], m)

return score + noise

def as_test_score(X):
# 点数を整数に
Y = X.round()
# 点数を0点から100点の間に制限
Y[Y <= 0] = 0
Y[Y >= 100] = 100
return Y

# 真の指標の数学,英語の比率
a = np.array([[1, 0.3]])
# 真の指標の中心
mu = np.array([56])
# 真の指標の分散
sigma = np.array([15])
# ノイズの中心を調整すると分布が上下左右にシフトする
noise_mu = np.array([0, 55])
# ノイズの分散を大きくすると数学と英語の点数の相関がなくなっていく
noise_sigma = np.array([5, 5])

# 1000人くらいほしい
X = as_test_score(make_dataset(a, mu, sigma, noise_mu, noise_sigma, 1000, 2))
# 平均点
m = np.mean(X, axis=0)

plt.figure(figsize=(16,8))
# ヒストグラムをプロット
plt.subplot(1,2,1)
plt.hist(X, bins=50, label=['Math', 'English'])
plt.legend()
plt.xlim([0,100])
# 散布図をプロット
plt.subplot(1,2,2)
plt.scatter(X[:,0], X[:,1], color='turquoise')
plt.quiver(0, 0, m[0], m[1], angles='xy', scale_units='xy', scale=1, color='purple')
plt.quiver(m[0], m[1], 20 * a[0][0], 20 * a[0][1], angles='xy', scale_units='xy', scale=1, color='blue')
plt.xlim([0,100])
plt.ylim([0,100])
plt.xlabel('Math')
plt.ylabel('English')
plt.show()


さて、ここで普通に学校で使われている「教科ごとの単純な合計点で生徒の優秀さを評価する」ことのどこに問題があるかを見ておきましょう。次元圧縮にもいろいろなやり方があるわけで、単純な合計点ももちろんひとつの指標となりえます。


正射影ベクトル

主成分分析と対比しながら単純な合計点の意味を理解するには、正射影(直交射影、Orthogonal Projection)の概念と結びつけて考えるとよいでしょう。正射影とは、直感的には「ある与えられた点についてあるベクトルに沿った情報を取り出す操作」です。

教科の単純な合計点$s$は、数学の点数を$x_1$,英語の点数を$x_2$とすると次の式で表すことができます。

$$

s = x_1 + x_2

$$これは、${\bf a}=(1,1)^{\mathrm T}, {\bf x}=(x_1, x_2)^{\mathrm T}$とおくと、

$$

s = {\bf a}^{\mathrm T} {\bf x} = 1\cdot x_1 + 1 \cdot x_2

$$というベクトルの内積として表現することができます。ここで観測ベクトル${\bf a}$との内積を取るとはどういうことか考察しておきます。いまもしも${\bf a}$が第一軸、すなわち${\bf a}={\bf e}_ 1=(1,0)^{\mathrm T}$であったならば、

$$

s = {\bf a}^{\mathrm T}{\bf x} = 1 \cdot x_1 + 0 \cdot x_2 = x_1

$$であり、これは数学の点数をそのまま取り出す操作です。もし${\bf a}={\bf e}_ 2=(0,1)^{\mathrm T}$であったならば英語の点数をそのまま取り出す操作になったでしょう。次の図で確認できる通り、$\|{\bf a}\| = 1$の観測ベクトル${\bf a}$との内積を取るということは、${\bf a}$に直交する光線(オレンジのベクトル)で点${\bf x}$の影を${\bf a}$に向けて写しとる操作になっています。

arrows.png

これを、${\bf x}$の${\bf a}$への正射影(あるいは直交射影)と言います。任意のベクトル${\bf x}$の任意の観測ベクトル${\bf a}$への正射影ベクトル${\bf p}$は以下の式で定義されます(『高校数学の美しい物語』にも解説があります)。

$$

{\bf p} = \frac{{\bf a}^{\mathrm T}{\bf x}}{\|{\bf a}\|^2}{\bf a}

$$定義は一見いかついですが、少し式を変形すると見通しがよくなります。

$$

\begin{align}

{\bf p} &= \frac{{\bf a}^{\mathrm T}{\bf x}}{\|{\bf a}\|^2}{\bf a} \\

&= \left(\frac{{\bf a}}{\|{\bf a}\|}\right)^{\mathrm T} {\bf x}\left(\frac{{\bf a}}{\|{\bf a}\|}\right) \\

&= ({\bf e}^{\mathrm T}{\bf x}){\bf e} \hspace{30px} where. {\bf e} = \left(\frac{{\bf a}}{\|{\bf a}\|}\right) \\

&= (\|{\bf e}\|\|{\bf x}\|cos\theta){\bf e} \\

& = (\|{\bf x}\|cos\theta){\bf e}

\end{align}

$$ここで$\theta$はベクトル${\bf a}$と${\bf x}$が成す角であり、${\bf e}$は${\bf a}$を${\bf a}$自身のノルムで割ったものですから${\bf a}$の方向を表す単位ベクトルです。この式の形ならば「${\bf a}$に直交する光線で点${\bf x}$の影を${\bf a}$に向けて写しとる操作」をイメージできるのではないでしょうか。

orthopro.png

この説明だけで正射影が「ある与えられた点についてあるベクトルに沿った情報を取り出す操作」であることをイメージするには十分かと思いますが、また少し見方を変えた説明もしておきましょう。${\bf x}$は一般には多次元のデータであり、たくさんの種類の情報を含んでいます。たとえば${\bf x}$が3次元のデータならば「${\bf e} _1, {\bf e} _2, {\bf e} _3$という3種類の(互いに直交するノルム$1$の)情報を、それぞれ$a _1, a _2, a _3$という量だけ含んでいる」ということであり、${\bf x}$を観測したときはそれらの情報の足し合わせとして観測されます(そう仮定してデータを扱います)。このときデータ${\bf x}$を数式で表現すると次のようになります。

$$

{\bf x} = a _1 {\bf e} _1 + a _2 {\bf e} _2 + a _3 {\bf e} _3

$$ここで${\bf x}$に対して${\bf e} _1^{\mathrm T}$との内積を取ってみましょう。内積は和に対して分配できるので、次のように計算できます。

$$

\begin{align}

{\bf e} _1^{\mathrm T} {\bf x} &= a _1 {\bf e} _1^{\mathrm T}{\bf e} _1 + a _2 {\bf e} _1^{\mathrm T}{\bf e} _2 + a _3 {\bf e} _1^{\mathrm T}{\bf e} _3 \\

&= a _1 \cdot 1 + a _2 \cdot 0 + a _3 \cdot 0 \\

&= a _1

\end{align}

$$つまり正射影を取ることで${\bf e} _1$と直交した(${\bf e} _1$とは関係のない)情報である${\bf e} _2, {\bf e} _3$を「削り落とす」ことができるのです。「ある与えられた点についてあるベクトルに沿った情報を取り出す操作」のイメージがつかめたでしょうか。本題からは少し逸れますが、ここで$a _1 = {\bf e} _1^{\mathrm T}{\bf x}$であることがわかり、同様に$a _2 = {\bf e} _2^{\mathrm T}{\bf x}, \, a _3 = {\bf e} _3^{\mathrm T}{\bf x}$であることが言えますから、

$$

{\bf x} = ({\bf e} _1^{\mathrm T}{\bf x}){\bf e} _1 + ({\bf e} _2^{\mathrm T}{\bf x}){\bf e} _2 + ({\bf e} _3^{\mathrm T}{\bf x}){\bf e} _3

$$と表すことができます。多次元に一般化すると${\bf x}$は

$$

{\bf x} = \sum _j ({\bf e} _j^{\mathrm T}{\bf x}){\bf e} _j

$$の形で表すことができ、これを${\bf x}$のフーリエ展開といいます。フーリエ展開は現代の信号処理やカーネル法の理論的な裏付けとなる関数解析という分野の重要な式のひとつであり「ベクトルに沿って情報を取り出す」という意味合いに変わりはありません。

さて、話を戻すことにしますと、各教科の単純な合計点を計算するための観測ベクトル${\bf a}$は、${\bf a}=(1,1)^{\mathrm T} = \sqrt{2}(\frac{1}{\sqrt{2}},\frac{1}{\sqrt{2}})^{\mathrm T}=\sqrt{2}{\bf e}$なので、単純な合計点$s$はベクトル${\bf e}=(\frac{1}{\sqrt{2}},\frac{1}{\sqrt{2}})^{\mathrm T}$に沿った正射影ベクトルの係数を$\sqrt{2}$倍したものになっていることがわかります。

$$

s = {\bf a}^{\mathrm T} {\bf x} = \sqrt{2}({\bf e}^{\mathrm T}{\bf x})=\sqrt{2}(\|{\bf x}\|cos\theta)

$$ここで$\theta$は${\bf e}$と${\bf x}$の成す角です。つまり単純な合計点で成績の優劣を比較することは「ベクトル${\bf e}=(\frac{1}{\sqrt{2}},\frac{1}{\sqrt{2}})^{\mathrm T}$への正射影ベクトルのノルムを比較する」こととたかだかスカラー倍の差を除いて等価なことがわかりました。

私たちの興味は「教科の単純な合計点で生徒の優秀さを評価する」ことのどこに問題があるか、ということでしたが、これはいま「得点ベクトル${\bf x}$をベクトル${\bf e}=(\frac{1}{\sqrt{2}},\frac{1}{\sqrt{2}})^{\mathrm T}$への正射影ベクトルのノルムで生徒の優秀さを評価する」ことのどこに問題があるかという問いに置き換えることができました。この指標は素朴で分かりやすいものですが、「数学の1点と英語の1点は同じ価値がある」という強い仮定を置いています。本当にそれでいいんでしょうか。私たちは


  • 数学の1点は英語の何点分に相当するか

ということを考えなくてよいのでしょうか?


単純な合計点で損をする人たち

では単純な合計点と主成分分析による評価で結果が食い違う例を確認しておきましょう。主成分分析を行って得られる観測ベクトルは$(\frac{1}{\sqrt{2}},\frac{1}{\sqrt{2}})^{\mathrm T}$に限らない、より一般の単位ベクトル${\bf e}$ですから、${\bf e} \neq (\frac{1}{\sqrt{2}},\frac{1}{\sqrt{2}})^{\mathrm T}$を適当に選んで${\bf e}_0 = (\frac{1}{\sqrt{2}},\frac{1}{\sqrt{2}})^{\mathrm T}$への正射影と結果が異なるケースが存在することを示せばよいです。

今回の主成分ベクトルは$(1, 0.3)^{\mathrm T}$なので、これに対応する単位ベクトルは${\bf e}=(\frac{1}{\sqrt{1.09}},\frac{0.3}{\sqrt{1.09}})^{\mathrm T}$です。「数学の$1$点は本当は英語の$\frac{1}{0.3}\approx 3.33$点に相当しますよ」というケースですね。「本来の${\bf e}$による評価が${\bf e}_ 0$による評価よりも小さくなってしまう範囲」を求めれば「損をする人たちの存在する範囲」が求まることになるので、次の不等式を解いてみましょう。

$$

{\bf e}^{\mathrm T}{\bf x} \lt {\bf e}_0^{\mathrm T}{\bf x}

$$$$

\therefore x _ 2 \lt \frac{\frac{1}{\sqrt{2}} - \frac{1}{\sqrt{1.09}}}{\frac{0.3}{\sqrt{1.09}} - \frac{1}{\sqrt{2}}} x _1 \approx 0.597 x _1

$$つまり下の図に引いた点線よりも下にいる人は本当の頭のよさよりもいささか不当な評価を受け、点線よりも上にいる人はラッキー。得してるわけです。実際に点線よりも上に行くほど得をして、点線よりも下に行くほど損をするため、うまく選ぶと下の図のように順位が入れ替わります。

[ここに図]

要するに「テストの問題の難易度が分野で異なる」という理由だけで不当な評価を受ける人がいるというわけです。

この節をこんな締め方をしていいものかわかりませんが、学業において重要なのは究極的には自分でどれだけ各分野を理解できたかです。自分のテストの点数を見て自分がどれだけその分野を理解できているかの指標にすることはある程度意味があるかもしれません。しかし総合順位はテストの製作者の影響で揺らいでしまうものであり、必ずしもあなたの努力や能力が公平に反映されるものではないのです。総合順位の上がり下がりに一喜一憂することなく、自分のペースで楽しく勉強ができるといいですね。


主成分分析の解き方

ここまで来ればあと理解すべきは数式だけです。主成分分析は「変換後の各確率変数の分散が最大になるような直交射影を求めよ」という問題です。すなわち

$$

z = {\bf a}^{\mathrm T}{\bf x} = a_1 x_1 + a_2 x_2

$$という新たな確率変数$z$を考えて、$z$の分散が最大となるような$a_1,a_2$を探せという最適化問題です。ただし制約条件として${\bf a}$は単位ベクトルですので、$a_1^2+a_2^2=1$を常に満たします(実はより一般に$a_1^2+a_2^2=r$としても以降の議論には影響を与えません)。

$z$の分散$Var(z)$はデータ数を$N$と置くと次のようになります。ちょっと計算が長いですが頑張ってください。ただし$K$は求めたいのが標本分散のとき$K=N$、母分散の不偏推定量のとき$K=N-1$となる定数です。

$$

\begin{align}

Var(z) &= \frac{1}{K} \sum_{i=1}^N (z_i -\overline{z})^2 \hspace{30px}where.\overline{z} = \frac{1}{N} \sum_{i=1}^N z_i\\

&= \frac{1}{K} \sum_{i=1}^N(z_i^2 - 2 z_i \overline{z} + \overline{z}^2) \\

&= \frac{1}{K} \sum_{i=1}^N\left\{(a_1 x_{1i} + a_2 x_{2i})^2 - 2 (a_1 x_{1i} + a_2 x_{2i}) (a_1 \overline{x}_1 + a_2 \overline{x}_2) + (a_1 \overline{x} _1 + a _2 \overline{x} _2) \right\} \\

&= \frac{1}{K} \sum _{i=1}^N\left\{a _1^2(x _{1i}^2 - 2 x _{1i}\overline{x} _1 + \overline{x} _1^2) + a _2^2(x _{2i}^2 - 2 x _{2i}\overline{x} _2 + \overline{x} _2^2) -2a _1 a _2 (x _{1i} - \overline{x} _1)(x _{2i} - \overline{x}) _2 \right\} \\

&= Var(x _1) a _1^2 + Var(x _2) a _2^2 + 2 Cov(x _1, x _2) a _1 a _2 \tag{1}

\end{align}

$$したがって$Var(z)$は$Var(x _1), Var(x _2), Cov(x _1, x _2)$という「計測されたデータから計算できる分散と共分散」を用いた$a _1, a _2$の2次式になっていることがわかります。このまま$Var(z)$を最大化しようとすると、$a _1, a _2$に正の大きな値を入れていけばいくらでも大きくできてしまいますから、$a _1, a _2$が取りうる値の範囲を$\|{\bf a}\|=1$に束縛します。

この等式制約付き最適化問題はラグランジュの未定乗数法によって解くことができます。ラグランジュの未定乗数法では目的関数$f(x)$を束縛条件$g(x)=0$のもとで最大化したいときにラグランジュ乗数$\lambda(\lambda > 0)$を用いて次のラグランジュ関数${\mathscr L}$(これはLの花文字です)に変換します。

$$

{\mathscr L}(x, \lambda) = f(x) - \lambda g(x)

$$この関数に所定の操作を行うことで最適解の候補を得ることができます。今回の記事はラグランジュの未定乗数法の解説記事ではありませんので詳しくは解説しませんが、ラグランジュの未定乗数法はよく登場する数理最適化の基本的な解法のひとつなので勉強しておくとよいでしょう。使い方はこれから説明いたします。

まず今回解くべき最適化問題を書き下します。

$$

\mathop{\rm arg\,max}\limits _{{\bf a}}Var(z)\hspace{30px}s.t. \|{\bf a}\|=1

$$数式に慣れていない人に説明するとこの数式は「$\|{\bf a}\|=1$という条件下で$Var(z)$を最大化する${\bf a}$を求めよ」という意味です($s.t.$は「such that(〜であるような)」の略ですが、数理最適化に関わる人は 「subject to(〜という条件下で)」と読むことを好みます。どちらでも同じ意味です)。このときラグランジュ関数は

$$

{\mathscr L}({\bf a}, \lambda) = Var(z) - \lambda (\|{\bf a}\| - 1)

$$となります(束縛条件は$\|{\bf a}\| = 1$のままでは使えないので$g(x)=0$の形にするため$g({\bf a})=\|{\bf a}\|-1 = 0$に変換します)。上の式はより具体的に書くと

$$

{\mathscr L}(a _1, a _2, \lambda) = Var(x _1) a _1^2 + Var(x _2) a _2^2 + 2 Cov(x _1, x _2) a _1 a _2 - \lambda (a _1^2 + a _2 ^2 - 1)

$$となります。ここから最適解の候補を得るための「所定の操作」をしていきます。やり方は簡単で各変数についての偏微分が$0$となる場所を探すだけです。すなわち

$$

\begin{cases}

\frac{\partial {\mathscr L}(a _1, a _2, \lambda)}{\partial a _1} &= 2Var(x _1) a _1+ 2Cov(x _1, x _2) a _2-2 \lambda a _1 = 0\\

\frac{\partial {\mathscr L}(a _1, a _2, \lambda)}{\partial a _2} &= 2Var(x _2) a _2+ 2Cov(x _1, x _2) a _1-2 \lambda a _2 = 0\\

\frac{\partial {\mathscr L}(a _1, a _2, \lambda)}{\partial \lambda} &= a _1^2 + a _2^2 - 1 = 0

\end{cases} \tag{1}

$$を満たすように連立方程式を解いて$a _1, a _2, \lambda$を定めます。ここで上式の第一式、第二式は行列形式で次のように書くことができます。

$$

\left(

\begin{matrix}

Var(x _1) & Cov(x _1, x _2) \\

Cov(x _1, x _2) & Var(x _2)

\end{matrix}

\right)\left(\begin{matrix}

a _1 \\

a _2

\end{matrix}

\right) = \lambda \left(\begin{matrix}

a _1 \\

a _2

\end{matrix}

\right) \tag{2}

$$ここで

$$

S = \left(

\begin{matrix}

Var(x _1) & Cov(x _1, x _2) \\

Cov(x _2, x _1) & Var(x _2)

\end{matrix}

\right) \tag{3}

$$とおくと$S$は確率変数$x _1, x _2$の分散共分散行列なので、$(2)$式は

$$

S {\bf a} = \lambda {\bf a} \tag{4}

$$を満たす${\bf a}$と$\lambda$を求めよ、という問題に変わります。$(1)$式の第三式と合わせれば、これは明らかに分散共分散行列$S$の固有値$\lambda$とノルムが$1$の固有ベクトル${\bf a}$を求めよという問題です(固有ベクトルのノルムには不定性があり、最初に束縛条件を一般に$\|{\bf a}\|=r$としておくと、ここでノルムが$r$の固有ベクトルを求めることになります)。

分散共分散行列は正定値対称行列なので$\lambda \gt 0$の固有値が2つと、それに対応する固有ベクトルが2本求まります(厳密には半正定値対称行列なので$\lambda = 0$となることもあります。もちろんそれは主成分の方向の分散が$0$である場合で、たとえばデータが1点に集まっていたり一直線に並んでいたりするとそうなります)。これが求まった「解の候補」であり、この2つのうち大きいほうの固有値に対応するのが第一主成分、小さいほうの固有値に対応するのが第二主成分です。第一主成分は最大化問題の大域最適解、第二主成分は最大化問題の局所最適解(極値)です。


主成分分析を使ってみる

そろそろ上にスクロールして問題を再確認するのが面倒でしょうから、同じ図をもう一度貼り付けておきます。

test_scores.png

私たちはこのようなデータから主成分ベクトル${\bf a}$(右図の青いベクトル)を求めたいのでした。では$(3),(4)$式あたりを使って主成分分析を実装していきましょう。まずは$(3)$式を用いてデータセットの分散共分散行列を計算します。


input

# rowvar は行方向、列方向のどちら向きに確率変数が格納されているかを指定する引数

S = np.cov(X, rowvar=False)
print(S)


output

[[ 258.67579479   69.76924625]

[ 69.76924625 43.78417518]]

numpy.covはデフォルトで母分散の不偏推定量を計算することに注意してください。すなわち$(1)$式で$K=N-1$の場合になっています。あとはこの固有値を求めればOKです。


input

lam, V = np.linalg.eig(S)

print(lam)
print(V)


output

[ 279.3406784    23.11929157]

[[ 0.95882613 -0.28399376]
[ 0.28399376 0.95882613]]

$V$の1列目${\bf v}_1 \fallingdotseq (0.96, 0.28)^{\mathrm T}$が第一主成分です。もともと数学$\colon$英語$=1\colon 0.3$としていますから計算結果とほぼ一致しています。これを図に描画して比較すると次のようになります(左図が真の${\bf a},$右図が推定された$\hat{\bf a}$)。

compare.png

あっさりできましたね。実際、主成分分析の計算自体は行列計算のライブラリさえあれば非常に簡単です。sklearnとか使えばもっと簡単ですが、それはまたの機会にしましょう。


N次元に一般化した場合

今までは入力$X$が2次元の場合について解説してきましたが、$X$が$N$次元のデータを持つ場合もまったく同じように計算することができます。$Cov(x _i, x _i) = Var(x _i)$とすれば、${\bf x}$が$N$次元のときの分散共分散行列$S$は$N \times N$の正方行列で、その$i$行$j$列目の成分は

$$

S _{ij} = Cov(x _j, x _i)

$$と表すことができます。あとは2次元のときと同様に

$$

S {\bf a} = \lambda {\bf a}

$$なる固有ベクトル${\bf a}$を求めるだけです。

せっかくなので主成分分析による次元圧縮の例を実装して今回の話を締めくくりましょう。

【注:追記予定】


終わりに

【注:追記予定】