前節ではscikit-learnを利用して主成分分析の考え方を学びました。ここでは、敢えてscikit-learnなしで主成分分析の計算のしくみを考えてみます。

同じ例題です。数学・理科・社会・英語・国語という5教科のテストの得点が1クラス20名分あったとして、5教科＝5次元を圧縮して、もっと少ない次元で端的に生徒の学力を評価するということをしてみたい。

5教科＝5つの変数を$x_{1}, x_{2}, x_{3}, x_{4}, x_{5}$として、求める主成分$z_{1}$の形を示します。

z = w_{1}x_{1} + w_{2}x_{2} + w_{3}x_{3} + w_{4}x_{4} + w_{5}x_{5}

各変数に係数$w_{1}, w_{2}, w_{3}, w_{4}, w_{5}$をつけた各成分の値を合計したものが主成分得点$z$となります。知りたいのは、この係数$w_{1}, w_{2}, w_{3}, w_{4}, w_{5}$をどういうふうな比率で割りふると情報量=分散が最大になるかということです。そこで、次のような条件をつけます。

{w_{1}}^2 + {w_{2}}^2 + {w_{3}}^2 + {w_{4}}^2 + {w_{5}}^2 = 1

係数をそれぞれ2乗した値を合計すると$1$になる、という条件のもとで情報量＝分散が最大となる配分比を探します。
このように、ある拘束条件のもとで多変数関数の最大値や最小値を求める方法としてラグランジュの未定乗数法があります。

ラグランジュの未定乗数法

変数$x, y$が拘束条件$g(x, y)=0$を満たしながら動くとき、関数$z=f(x, y)$が最大値･最小値となる$x, y$には次の式が成り立つ。

F(x, y, λ) = f(x, y) - λg(x, y)　とすると、

\displaystyle \frac{∂F}{∂x} = 0,　\frac{∂F}{∂y} = 0,　\frac{∂F}{∂λ} = 0

以上が最もシンプルな2変数、1条件でのラグランジュの未定乗数法の考え方です。$λ$（ラムダ）という記号ですが、ここでは「これは定数です」という理解にとどめます。また$∂$という記号は、ラウンドとかパーシャルとかデルとかいくつも呼び名がありますが、偏微分を表します。
一つ例に引けば、$\displaystyle \frac{∂F}{∂x} = 0$ は、関数$F$を変数$x$で偏微分した値が$0$という意味です。

偏微分

まず、ふつうの微分の意味を確認します。
$x, y$平面上に$y=f(x)$という関数の曲線を描いてみます。
$x$座標上の$x$という点に注目して、この点での高さを調べると、$y=f(x)$上の点だから$f(x)$です。また$x$から少しずれた$x+Δx$という点に注目してみます。この点での高さは関数$y=f(x)$に$x+Δx$を当てはめたものだから$y$座標は$f(x+Δx)$となります。
そして点と点との傾きを出します。傾きは$x$の変化量分の$y$の変化量です。分母の$x$の変化量は$Δx$、分子の$y$の変化量は後引く前で$f(x+Δx)-f(x)$となります。

微分とは何かというと、この$Δx$を$0$に近づけていくことでした。つまり$x+Δx$が極限まで$x$に近づくと、$x$という点でのその一瞬の傾きを表わす接線が引けます。この傾きが微分の意味です。

\frac{dy}{dx}=\lim_{\Delta x \to 0} \frac{f(x+\Delta x)-f(x)}{\Delta x}

$x$が少し変化すると$y$はどれだけ変化するか、微分とは、$x$の微小変化に対する$y$の変化の割合です。

ここまでは、$x$が1つ決まれば$y$が決まる、1変数関数の場合の微分でした。しかし2つの値が決まって初めてその値が決まる2変数関数$f(x,y)$もあれば、3つの値が決まって初めてその値が決まる3変数関数$f(x,y,z)$もあります。一般には多変数関数というものがあって、実は多変数関数の微分を考えるのが偏微分です。
2変数関数のグラフを描いてみます。

$x$と$y$が決まると、高さ$z$が決まります。
下の平面上のどの点にも対応する高さ$z=f(x,y)$があって、曲面のような感じになります。そして、$x$軸方向にちょっと進むのか$y$軸方向にちょっと進むのかで高さ$z$は変わってきます。ちょっと動いたときに高さがどう変わるのかは、大きく2つの方向があります。それが$x$軸ではどうなのか、$y$軸ではどうなのかというのが偏微分です。

$y$を固定して$x$がちょっとだけ動いたときの傾き、つまり$x$軸方向の接線を引いたものが青色の線です。式に表わします。

\frac{∂z}{∂x}=\lim_{\Delta x \to 0} \frac{f(x+\Delta x, y)-f(x, y)}{\Delta x}

$y$は固定されています。$x$方向にちょっとだけ動いたときに高さ$z$がどれだけ動くのか、$x$の微小変化に対する$z$の変化の割合が偏微分です。
したがって$y$の偏微分は、

\frac{∂z}{∂y}=\lim_{\Delta y \to 0} \frac{f(x, y + \Delta y)-f(x, y)}{\Delta y}

$x$が固定されて、$y$方向にちょっとだけ動いたときに高さ$z$がどれだけ動くのか、$y$の微小変化に対する$z$の変化の割合で、赤い線がそれです。

さて、偏微分した値が$0$ということの意味を考えてみます。
関数$z=f(x,y)$で、変数$x$をちょっとだけ動かすと$z$の変化が$0$である。つまり変化がない、それ以上大きくも小さくもならないとすれば、そこが最大か最小だからです。
さらに、ラグランジュの未定乗数法が示すように、同時に$x$も$y$も偏微分した値が$0$であるとき、2変数関数の曲面上で真の山頂（最大値）か真の谷底（最小値）だということになります。

偏微分の方程式から分散共分散行列を導く

5教科の例題に戻ります。
くり返しになりますが、主成分のかたちは次のとおりです。
$p=ax+by+cu+dv+ew$

主成分$p$の係数$a,b,c,d,e$は定数で、この$p$の分散が最大になるように決められます。なお、知りたいのは定数$a,b,c,d,e$の比率のバランスなので$a,b,c,d,e$には次の条件がつきます。
$a^2+b^2+c^2+d^2+e^2=1$

最大値を求める主成分の分散${s_{p}}^2$を定義します。個体数を$n$として次のように表わされます。
$\displaystyle Sp^2=\frac{1}{n} \{ {(p_{1}-\bar{p})}^2+{(p_{2}-\bar{p})}^2+...+{(p_{n}-\bar{p})}^2\}$

ここで、ラグランジュの未定乗数法を適用し、関数$L$を定義します。
$\displaystyle L=\frac{1}{n} \{ {(p_{1}-\bar{p})}^2 + {(p_{2}-\bar{p})}^2 + ... + {(p_{n}-\bar{p})}^2 \} - λ(a^2 + b^2 + ... + e^2 - 1)$
とおくと、次の方程式が成り立ちます。
$\displaystyle \frac{∂L}{∂a}=0 , \frac{∂L}{∂b}=0 , \frac{∂L}{∂c}=0 , \frac{∂L}{∂d}=0 , \frac{∂L}{∂e}=0$

試みに一番目の式に注目します。
$\displaystyle \frac{∂L}{∂a}=0$

関数$L$の式を分子にもってきて、分母に$a$をおいて偏微分します。$a$以外は固定といいましたが、具体的には定数とみなして$b,c,d,e$は消してしまいます。
$\displaystyle \frac{∂L}{∂a} = \frac{2}{n} \{ (p_{1}-\bar{p})(x_{1}-\bar{x}) + (p_{2}-\bar{p})(x_{2}-\bar{x}) + ... + (p_{n}-\bar{p})(x_{n}-\bar{x}) \} - 2λa = 0$

主成分$p$の式をもってきて｛｝の中を展開し、分散と共分散の定義をつかって変形します。
$\displaystyle 2 \{ a{s_{x}}^2 + bs_{xy} + cs_{xu} + ds_{xv} + es_{xw} \} - 2λa = 0$

さらに整理すると次の式が得られます。
$a{s_{x}}^2 + bs_{xy} + cs_{xu} + ds_{xv} + es_{xw} = λa$

同様の式が残りの$b,c,d,e$からも得られます。
$as_{xy} + b{s_{y}}^2 + cs_{yu} + ds_{yv} + es_{yw} = λb$
$as_{xu} + bs_{yu} + c{s_{u}}^2 + ds_{uv} + es_{uw} = λc$
$as_{xv} + bs_{yv} + cs_{uv} + d{s_{v}}^2 + es_{vw} = λd$
$as_{xw} + bs_{yw} + cs_{uw} + ds_{vw} + e{s_{w}}^2 = λe$

この5つの式を行列にまとめて整理すると次のようになります。

\begin{pmatrix}
{Sx}^2 & S_{xy} & S_{xu} & S_{xv} & S_{xw}\\
S_{xy} & {Sy}^2 & S_{yu} & S_{yv} & S_{yw}\\
S_{xu} & S_{yu} & {Su}^2 & S_{uv} & S_{uw}\\
S_{xv} & S_{yv} & S_{uv} & {Sv}^2 & S_{vw}\\
S_{xw} & S_{yw} & S_{uw} & S_{vw} & {Sw}^2
\end{pmatrix}
\begin{pmatrix}
a\\
b\\
c\\
d\\
e
\end{pmatrix}
=
λ
\begin{pmatrix}
a\\
b\\
c\\
d\\
e
\end{pmatrix}

この左辺の正方行列（行と列の要素数が同じ行列）は分散共分散行列になっています。この式を満たす列ベクトル$a,b,c,d,e$を固有ベクトル、$λ$を固有値と呼びます。この固有値$λ$は主成分の分散になっています。

⑴ ライブラリをインポートする

import pandas as pd
import numpy as np
import matplotlib.pyplot as plt

⑵ データを作成する

# 生徒20名の5教科の得点を作成
arr = np.array([[71,64,83,100,71], [34,48,67,57,68], [58,59,78,87,66], [41,51,70,60,72],
                [69,56,74,81,66], [64,65,82,100,71], [16,45,63,7,59], [59,59,78,59,62],
                [57,54,84,73,72], [46,54,71,43,62], [23,49,64,33,70], [39,48,71,29,66],
                [46,55,68,42,61], [52,56,82,67,60], [39,53,78,52,72], [23,43,63,35,59],
                [37,45,67,39,70], [52,51,74,65,69], [63,56,79,91,70], [39,49,73,64,60]])

# カラム名を指定してデータフレームを作成
df = pd.DataFrame(data = arr,
                  columns = ['数学', '理科', '社会', '英語', '国語'])

⑶ 分散共分散行列を求める

# データフレームをNumpyのarray型に変換
arr = np.array(df)

# 分散共分散行列を求める
cov = np.cov(arr, rowvar=0, bias=0)

分散共分散行列（variance covariance matrix）は、Numpyのcov関数をつかって求めます。
numpy.cov(data, rowvar=0, bias=0)のかたちで、引数のbias=0は不偏分散を指定しています（標本分散ならbias=1）。rowvar=0はデータが縦方向に並んでいる場合で、横方向の場合はrowvar=1を指定します。

⑷ 固有値･固有ベクトルを求める

# 固有値･固有ベクトルを求める
eig_val, eig_vec = np.linalg.eig(cov)

print("固有値：\n", eig_val) 
print("固有ベクトル：\n", eig_vec)

Numpyのlinalg.eig()関数は、行列を指定するだけで固有値と固有ベクトルの2つの配列をタプルで返します。タプルとは、複数のデータがまとめられたデータ型で、いったん生成した後に変更ができないオブジェクトです。
ここでは、固有値をeig_val、固有ベクトルをeig_vecとして別々に格納しました。

固有値のデータの並びに注意してください。scikit-learnで見たような寄与率の順には並んでいません。そこで並び替えをして、データフレームに変換して見やすくしていきます。

# 固有値･固有ベクトルを1つのarrayに連結
eigen = np.vstack((eig_val,eig_vec))

# 先頭行を基準にソート
eigen[:, eigen[0,:].argsort()[::-1]]

Numpyのvstack((a,b))関数で縦方向（axis=0）に連結します。
Numpyのargsort関数で特定の行を基準にソートします。このとき[::-1]]をつけて降順を指定しています。

# 固有値をデータフレームに変換
eval_df = pd.DataFrame(data = eigen[0,:].reshape([1, 5]),
                       columns = ["主成分{}".format(x+1) for x in range(len(df.columns))],
                       index = ["固有値"])
eval_df

# 主成分負荷量をデータフレームに変換
evec_df = pd.DataFrame(data = eigen[1:,:].reshape([5, 5]),
                       columns = ["主成分{}".format(x+1) for x in range(len(df.columns))],
                       index = ['数学', '理科', '社会', '英語', '国語'])
evec_df

いわゆる固有値問題、線形代数で有名な問題に行き着くのです。ここで固有値問題について深めたいところですが、先へ進みます。次節は因子分析をscikit-learnを利用して行ないます。