この記事は、私が参加しているオンラインセミナー「統計検定準一級勉強会」の講義ノートです。セミナーで使用しているテキスト「統計学実践ワークブック」は統計学のエッセンスを簡潔にまとめている良書ですが、数式などは結果のみで導出方法が記載されていないので、この記事で、できるだけフォローしました。このテキストを使用して勉強している人の理解の手助けに少しでもなれば幸いです。間違いや助言があればご指摘していただくとありがたいです。
$$
\newcommand{\argmax}{\mathop{\rm arg~max}\limits}
\newcommand{\argmin}{\mathop{\rm arg~min}\limits}
$$
23 判別分析（その２）
正準判別分析（P.205）
正準判別分析は、重判別分析ともよばれ1)、サンプルが属する群が 2 より多い場合の多群の判別において効率的な判別の方法です。教科書のこの節の冒頭ではマハラノビス平方距離による多群の判別に触れていますが、参考文献 1にも書かれているように、
例えば5群を判別したいとすると、2群ごとの判別分析を10種類行う必要があります。
というように、相当複雑な作業を必要とします。そこで考案されたのがこの正準判別分析です。正準判別分析は、フィッシャーの判別分析を多群に拡張したようなもので、多群のデータを判別性が高くなるように低次元に射影する方法です。参考文献 1 では、これを次のように表現しています。
まず最初に多群の中で最も離れている2群のプロットの重なりが最小になるような判別軸1――図18.2.1では正準軸1――を見つけて、その軸に多群のプロットを正射影して判別スコア1を求めます。 次に、その軸では重なりが大きかった別の2群の重なりが最小になるような判別軸2を見つけて、今度はその軸に多群のプロットを正射影して判別スコア2を求めます。 こうして多群を複数の判別軸で判別すると、最終的には(群数−1)個の判別軸を見つけ出し、判別軸の数だけ判別スコアを計算することができます。
言葉だけの説明はわかりにくいので、是非、参考文献 1 の「図18.2.1 重判別分析の概念図」および「図18.2.2 正準スコアのプロット」を参照しながらこの文章を読んでください。
正準判別分析の定式化
では、さっそく教科書に沿って正準判別分析の定式化を行います。教科書は結果だけを書いていますが、ここではできるだけ計算式を詳細にフォローしていきます。
まず、群の数を $g$ とし、群 $G_j,\hspace{3pt} j=1,\cdots,g$ に属するサンプルを $\lbrace \boldsymbol{x}_i^{(j)} \rbrace_{i=1}^{n_j},\hspace{3pt} j=1,\cdots,g$ とします。群 $G_j$ の平均ベクトルを
$$
\bar{\boldsymbol{x}}^{(j)} = \frac{1}{n_j} \sum_{i=1}^{n_j} \boldsymbol{x}_i^{(j)}
\quad・・・\quad(1)
$$
とし、分散共分散行列を
$$
S_j = \frac{1}{n_j - 1} \sum_{i=1}^{n_j}
\left( \boldsymbol{x}_i^{(j)} - \bar{\boldsymbol{x}}^{(j)} \right)
\left( \boldsymbol{x}_i^{(j)} - \bar{\boldsymbol{x}}^{(j)} \right)^{\mathrm{T}}
\quad・・・\quad(2)
$$
とし、全サンプルの平均ベクトルを
$$
\bar{\boldsymbol{x}}
= \left( \sum_{j=1}^{g} n_j \right)^{-1}
\sum_{j=1}^g \sum_{i=1}^{n_j} \boldsymbol{x}_i^{(j)}
= \frac{1}{n_1 + \cdots + n_g}
\sum_{j=1}^g n_j \bar{\boldsymbol{x}}^{(j)}
\quad・・・\quad(3)
$$
とします。全サンプルの平均ベクトルは、フィッシャーの線形判別の場合と違って各群のサンプル数 $n_j$ の重み付き平均になっていることに注意してください。フィッシャーの線形判別では、P204 の下から3行目の式（判別分析（その１）の(25)式）のように、中点を$\frac{\left( \bar{\boldsymbol{x}}^{(1)} + \bar{\boldsymbol{x}}^{(2)} \right)}{2}$ としています。全サンプルの平均ベクトルであれば $\frac{\left( \bar{n_1 \boldsymbol{x}}^{(1)} + n_2 \bar{\boldsymbol{x}}^{(2)} \right)}{n_1 + n_2} $ とすべきです。なぜフィッシャーの線形判別では全サンプルの平均ベクトルではなく中点としているのかはわかりません。
射影ベクトルを求める
さて、次に射影ベクトル $\boldsymbol{w}$ を求めましょう。やり方はフィッシャーの線形判別の場合と同じで、群間変動を最大化し、かつ、群内変動を最小化する射影方向を求めます。ただし、後でわかるように、この $\boldsymbol{w}$ はたった1つではありません。$g$ 個の群を判別するためには、$g - 1$ 個の $\boldsymbol{w}$ が必要です。それはさておき、まず、群間変動を求めるために、群間変動行列 $S_B$ を次式で定義します。
$$
S_B = \sum_{j=1}^{g} n_j
\left( \bar{\boldsymbol{x}}^{(j)} - \bar{\boldsymbol{x}} \right)
\left( \bar{\boldsymbol{x}}^{(j)} - \bar{\boldsymbol{x}} \right)^{\mathrm{T}}
\quad・・・\quad(4)
$$
ここでも、フィッシャーの線形判別の場合と違います。フィッシャーの線形判別の場合は、教科書P204の上から6行目の式（判別分析（その１）の(6)式）にあるように、各群のサンプル数は考慮されていません。
次に、群内変動です。群内変動を表す群内変動行列 $S_W$ を次のように定義します。
$$
S_W = \sum_{j=1}^{g} \left( n_j - 1 \right) S_j
\quad・・・\quad(5)
$$
これは、教科書P204の上から11行目の式（判別分析（その１）の(11)式）と、係数 $\frac{1}{n_1 + n_2 -2}$ を除いて一致します。その理由としては、群内変動は、各群のサンプル $\boldsymbol{x}_i^{(j)}$ とその平均 $\bar{\boldsymbol{x}}^{(j)}$ の距離の二乗和なので、サンプル数は既に織り込み済みだからでしょう（(2)式の $\sum_{i=1}^{n_j}$ の $i$ の上限 $n_j$）。群間変動と群内変動のこの一貫性のなさは不思議です。
さて、ここまできたら、あとはフィッシャーの線形判別と同様に射影軸上での群間、群内の分散比 $\lambda({\boldsymbol{w}})$ を最大化する極値問題を解くだけです（判別分析（その１）の(4)式）。
射影軸上での群間分散は、(4)式の群間変動行列 $S_B$ を用いて次のように表されます。
$$
群間分散 = \frac{1}{g - 1} \boldsymbol{w}^{\mathrm{T}} S_B \boldsymbol{w}
\quad・・・\quad(6)
$$
これは、判別分析（その１）の(7)式の計算方法で証明できます。判別分析（その１）の(7)式は、群内分散を求めていますが、群間、群内に関係なくこのような流れで証明できます。注目していただきたいのは、結果の式 $\boldsymbol{w}^{\mathrm{T}} S \boldsymbol{w}$ です。分散共分散行列 $S$ を $\boldsymbol{w}^{\mathrm{T}}$ と $\boldsymbol{w}$ で挟むと $\boldsymbol{w}$ 方向に射影した分散が得られるということです。サンプルを $\boldsymbol{x}$ とすると、 $\boldsymbol{w}^{\mathrm{T}} \boldsymbol{x}$ は、サンプルを $\boldsymbol{x}$ の $\boldsymbol{w}$ 方向への射影であることを思い出すと、何か $S \boldsymbol{w}$ が"分散ベクトル"のような感じがして面白いです（ただ、どうしてそうなるのかはわかりませんが・・・）。
さて、同様にして射影軸上での群内分散を次式で求めることができます。
$$
群内分散 = \frac{1}{n_1 + \cdots + n_g - g} \boldsymbol{w}^{\mathrm{T}} S_W \boldsymbol{w}
\quad・・・\quad(7)
$$
これは、先ほど触れた教科書P204の上から11行目の式（判別分析（その１）の(11)式）と群数が $g$ 個あることを除いて一致します。
以上から、群間、群内の分散比 $\lambda{(\boldsymbol{w})}$ は次式になります。
$$
\lambda{(\boldsymbol{w})} = \frac{n_1 + \cdots + n_g - g}{g - 1} \frac{\boldsymbol{w}^{\mathrm{T}} S_B \boldsymbol{w}}{\boldsymbol{w}^{\mathrm{T}} S_W \boldsymbol{w}}
\quad・・・\quad(8)
$$
これに適当な制約条件を導入して、いつものようにラグランジュの未定乗数法で解けばが $\boldsymbol{w}$ が求まります。これについては、参考文献 1 に詳細な数式展開が書いてあるので読んでみてください。ここでは、別の方法で $\boldsymbol{w}$ を求めてみます。
(8)式を以下のように変形します。
$$
\boldsymbol{w}^{\mathrm{T}} S_B \boldsymbol{w} = \frac{g - 1}{n_1 + \cdots + n_g - g} \lambda{(\boldsymbol{w})} \boldsymbol{w}^{\mathrm{T}} S_W \boldsymbol{w}
\quad・・・\quad(9)
$$
ここで、あらためて
$$
\frac{g - 1}{n_1 + \cdots + n_g - g} \lambda{(\boldsymbol{w})} \stackrel{\mathrm{def}}{=} \lambda
\quad・・・\quad(10)
$$
とおくと、(9)式は次のようになります。
$$
\boldsymbol{w}^{\mathrm{T}} S_B \boldsymbol{w} = \lambda \boldsymbol{w}^{\mathrm{T}} S_W \boldsymbol{w}
\quad・・・\quad(11)
$$
両辺に左から $\boldsymbol{w}^{\mathrm{T}}$ を掛けているのを消去すると
$$
S_B \boldsymbol{w} = \lambda S_W \boldsymbol{w}
\quad・・・\quad(12)
$$
さらに、両辺に $S_W^{-1}$ を掛けると
$$
S_W^{-1}S_B \boldsymbol{w} = \lambda \boldsymbol{w}
\quad・・・\quad(13)
$$
となり、これは、$\boldsymbol{w}$ が行列 $S_W^{-1}S_B$ の固有値 $\lambda$ に対応する固有ベクトルになっていることを示しています。でも、不思議です。ラグランジュの未定乗数法で極値問題を解かなくても同じ結果が得られるとは。きっと、固有ベクトル $\boldsymbol{w}$ が(8)式の $\lambda(\boldsymbol{w})$ の極大値を与えるのでしょう。そして、そのときの $\lambda(\boldsymbol{w})$ が、(10)式を経て固有値になっているのでしょう。
固有ベクトルの数
教科書P.206の上から8行目から10行目にかけて次の記述がみられます。
行列 $S_W^{-1}S_B$ の $d$ 個の固有値、固有ベクトルを $(\lambda_k, \boldsymbol{w}_k), \hspace{3pt}k=1,\cdots,d$ （ただし $\lambda_1 \ge \lambda_2 \ge \cdots \lambda_d$）とする。 $S_B$ と $\bar{\boldsymbol{x}}$ の定義から行列 $S_B$ のランクはたかだか $g - 1$ であり、$d \le \min{\lbrace p, g - 1 \rbrace}$ であることに注意する。
これは、何を言っているのでしょう？「ランク」という用語が唐突に出てきました。どうやらランクは固有値の数に関係がありそうです。
参考資料 2 によれば、行列のランクとは「行列にデータがどれほど詰まっているか」ということだそうです。もう少し具体的に言えば、「行列のランクは、列の線形独立なベクトルの最大個数」ということです。線形独立なベクトルって何でしょう？線形独立という用語の反対語は線形従属です。いま、2 つのベクトル $\boldsymbol{a}$ と $\boldsymbol{b}$ があったとしましょう。任意のスカラーの定数 $c \ne 0$ に対して $\boldsymbol{a}$ と $\boldsymbol{b}$ の間に $\boldsymbol{a} = c \boldsymbol{b}$ の関係が成り立つとき、$\boldsymbol{a}$ と $\boldsymbol{b}$ は線形従属といいます。$\boldsymbol{a}$ が $\boldsymbol{b}$ を定数倍して得られるのであれば、$\boldsymbol{a}$ の情報量はないというわけです（$\boldsymbol{b}$ があれば十分ということでしょうか）。逆に、$\boldsymbol{a} = c \boldsymbol{b}$ という関係が成り立たない場合、$\boldsymbol{a}$ と $\boldsymbol{b}$ は線形独立といいます。このように、行列の列ベクトルのうち、線形独立なベクトルが何個あるか、つまり、この行列にどれだけの情報量があるかを表す指標がランクなわけです。
ランクには次の性質があります。
- $A$ を $m \times n$ 行列とすると、$\mathrm{rank}(A) \le \min{(m, n)} $
- $\mathrm{rank}(AB) \le \min{(\mathrm{rank}(A), \mathrm{rank}(B))} $
- $\mathrm{rank}(A + B) \le \mathrm{rank}(A) + \mathrm{rank}(B) $
- $\boldsymbol{0}$ でない任意のベクトルを $\boldsymbol{a}, \boldsymbol{b}$ として、ベクトル $\boldsymbol{a}$ と $\boldsymbol{b}$ の外積でできた行列 $A = \boldsymbol{a} \boldsymbol{b}^{\mathrm{T}}$ に対して $\mathrm{rank}(A) = 1$
では、この性質を利用して、「行列 $S_B$ のランクはたかだか $g - 1$ である」ことを示しましょう。(4)式より
$$
\begin{eqnarray}
\mathrm{rank} (S_B) &=& \mathrm{rank} \left( \sum_{j=1}^{g} n_j
\left( \bar{\boldsymbol{x}}^{(j)} - \bar{\boldsymbol{x}} \right)
\left( \bar{\boldsymbol{x}}^{(j)} - \bar{\boldsymbol{x}} \right)^{\mathrm{T}} \right) \\
&\le& \sum_{j=1}^{g} \mathrm{rank} \left( \left( \bar{\boldsymbol{x}}^{(j)} - \bar{\boldsymbol{x}} \right)
\left( \bar{\boldsymbol{x}}^{(j)} - \bar{\boldsymbol{x}} \right)^{\mathrm{T}} \right) \\
&\le& \sum_{j=1}^{g} 1 = g
\quad・・・\quad(14)
\end{eqnarray}
$$
となります。ここで、1行目から2行目に移るときランクの性質 3 を使い、2行目から3行目に移るときランクの性質 4 を使いました。
ところで、(14)式に現れる $\bar{\boldsymbol{x}}^{(j)} - \bar{\boldsymbol{x}}$ ですが、この中の $\bar{\boldsymbol{x}}$ は、(3)式で与えられる全サンプルの平均ベクトルです。中心化を行うために $\bar{\boldsymbol{x}}^{(j)}$ から $\bar{\boldsymbol{x}}$ を引いているのですが、そのために $\bar{\boldsymbol{x}}^{(j)} - \bar{\boldsymbol{x}}$ は、すべての $j$ について線形独立にはなっていません。たとえば、$\bar{\boldsymbol{x}}^{(g)} - \bar{\boldsymbol{x}}$ は $\bar{\boldsymbol{x}}^{(1)} \sim \bar{\boldsymbol{x}}^{(g - 1)}$ から次のように計算することができます。
$$
\bar{\boldsymbol{x}}^{(g)} - \bar{\boldsymbol{x}} = (n - 1) \bar{\boldsymbol{x}} - \sum_{j=1}^{g - 1} \bar{\boldsymbol{x}}^{(j)}
\quad・・・\quad(15)
$$
したがって、行列 $\left( \bar{\boldsymbol{x}}^{(g)} - \bar{\boldsymbol{x}} \right)
\left( \bar{\boldsymbol{x}}^{(g)} - \bar{\boldsymbol{x}} \right)^{\mathrm{T}}$ の情報量はゼロになります。よって、(14)式は $\mathrm{rank} (S_B) \le g - 1$ となり、「行列 $S_B$ のランクはたかだか $g - 1$ である」ことを示すことができました。
次に、行列 $S_W^{-1}S_B$ の固有ベクトルの数 $d$ が、$d \le \min{(p, g - 1)}$ という関係を満たすことを示します。
行列 $A$ のランクと $A$ の非ゼロの固有値の数の関係については、
$$
「A の（重複度も含めた）非ゼロの固有値の数」 \le \mathrm{rank} (A)
\quad・・・\quad(16)
$$
という関係があります3)。(16)式の $A$ に $S_W^{-1}S_B$ を代入すれば
$$
d \le \mathrm{rank} (S_W^{-1}S_B) \le \min{\left(\mathrm{rank} (S_W^{-1}), \mathrm{rank} (S_B) \right)}
\quad・・・\quad(17)
$$
となります。ここで、ランクの性質 2 を使いました。
まず、$\mathrm{rank} (S_W^{-1})$ について考えてみましょう。分散共分散行列 $S_W$ は正則なので、$\mathrm{rank} (S_W^{-1}) = \mathrm{rank} (S_W)$ です。(5)式から
$$
\mathrm{rank} (S_W) = \mathrm{rank} \left( \sum_{j=1}^{g} \left( n_j - 1 \right) S_j \right) \le \sum_{j=1}^{g} \mathrm{rank} (S_j)
\quad・・・\quad(18)
$$
となります。ここで、ランクの性質 3 を利用しました。また、(2)式から
$$
\begin{eqnarray}
\mathrm{rank} (S_j)
&=& \mathrm{rank} \left( \frac{1}{n_j - 1} \sum_{i=1}^{n_j}
\left( \boldsymbol{x}_i^{(j)} - \bar{\boldsymbol{x}}^{(j)} \right)
\left( \boldsymbol{x}_i^{(j)} - \bar{\boldsymbol{x}}^{(j)} \right)^{\mathrm{T}} \right) \\
&\le& \sum_{i=1}^{n_j} \mathrm{rank} \left( \left( \boldsymbol{x}_i^{(j)} - \bar{\boldsymbol{x}}^{(j)} \right)
\left( \boldsymbol{x}_i^{(j)} - \bar{\boldsymbol{x}}^{(j)} \right)^{\mathrm{T}} \right) \\
&\le& n_j - 1
\quad・・・\quad(19)
\end{eqnarray}
$$
となります。ここで、1行目から2行目に移るとき、ランクの性質 3 を利用し、2行目から3行目に移るとき、$\mathrm{rank} (S_B)$ の導出のプロセスと同様に行いました。(18)式と(19)式から
$$
\mathrm{rank} (S_W) \le \sum_{j=1}^{g} \mathrm{rank} (S_j) \le \sum_{j=1}^{g} (n_j - 1) = n - g
\quad・・・\quad(20)
$$
が得られます。ここで、$n$ は、全サンプル数 $\sum_{j=1}^{g} n_j$ です。
次に、$S_B$ は$p \times p$ の行列なので、ランクの性質 1 より、 $\mathrm{rank} (S_B) \le p$ です。したがって、 前述した $\mathrm{rank} (S_B) \le g - 1$ を考慮すると、
$$
\mathrm{rank} (S_B) \le \min{(p, g - 1)}
\quad・・・\quad(21)
$$
となります。
(17)式に(20)式と(21)式を代入して、$\mathrm{rank} (S_W^{-1}) = \mathrm{rank} (S_W)$ を考慮すると
$$
d \le \min{(n - g, \min{(p, g - 1)})}
$$
を得ます。ここで、全サンプル数 $n$ を十分にとって $n - g \ge \min{(p, g - 1)}$ とすれば、最終的に $d \le \min{(p, g - 1)}$ が成立します。これで、本節の冒頭で掲げた教科書の記述を証明することができました。
2次元平面にプロット
射影軸（固有ベクトル $\boldsymbol{w}_k$ 方向の直線）を求めることができたので、サンプルを射影軸に射影して可視化を試みましょう。ここでの作業は、主成分分析の主成分得点のプロットとほぼ同じです。「射影ベクトルを求める」の節でも述べたように、固有ベクトル $\boldsymbol{w}_k$ に対応する固有値 $\lambda_k$ は(8)式で示す群間、群内分散比 $\lambda (\boldsymbol{w})$ に比例するので、その分散比がもっとも大きい（つまり、固有値がもっとも大きい）固有ベクトル $\boldsymbol{w}_1$ と2番目に大きい固有ベクトル $\boldsymbol{w}_2$ の方向に座標軸をとり（それぞれ、横軸、縦軸）、その座標平面上に $\left( \boldsymbol{w}_1^{\mathrm{T}} \boldsymbol{x}_i, \boldsymbol{w}_2^{\mathrm{T}} \boldsymbol{x}_i \right), \hspace{3pt} i = 1, \cdots, n$ をプロットします。ここで、$\boldsymbol{x}_i$ は $i$ 番目のサンプルです。では、さっそくやってみましょう。
Step 1. データ読み込み
使用するデータは、参考資料 4 の「表18.1.1 3群の検査項目」です。これをコピペして、TSVファイルに加工します。といっても、エクセルにペーストして、それをテキストファイルにコピーして保存するだけです。ファイル名を "labo_data.tsv" とします。
# データ読み込み
# データの出典
# 18.1 多群の判別 の 表18.1.1 3群の検査項目
# http://www.snap-tck.com/room04/c01/stat/stat18/stat1801.html
import pandas as pd
labo_data = 'labo_data.tsv'
T_data = pd.read_csv(labo_data, sep='\t')
TSVファイルをデータフレーム T_data に読み込みます。下記に読み込んだ T_data を示します。サンプルは全部で 37 個あり（$n = 37$）、3つの群（$g = 3$）にわかれています。表の列「群」に群ラベル（「正常」, 「疾患A」, 「疾患B」）が記入されています。各々のサンプル数は、それぞれ、$n_1 = 10, n_2 = 15, n_3 = 12$です。また、サンプルデータの特徴量（要因）は、「検査項目1」～「検査項目5」の 5 つ（$p = 5$）です。
| ID | 群 | 検査項目1 | 検査項目2 | 検査項目3 | 検査項目4 | 検査項目5 | |
|---|---|---|---|---|---|---|---|
| 0 | N01 | 正常 | 0 | 2 | 2 | 4 | 3 |
| ・・・ | ・・・ | ・・・ | ・・・ | ・・・ | ・・・ | ・・・ | ・・・ |
| 9 | N10 | 正常 | 0 | 8 | 6 | 2 | 5 |
| 10 | A01 | 疾患A | 0 | 1 | 3 | 6 | 5 |
| ・・・ | ・・・ | ・・・ | ・・・ | ・・・ | ・・・ | ・・・ | ・・・ |
| 24 | A15 | 疾患A | 0 | 7 | 8 | 5 | 5 |
| 25 | B01 | 疾患B | 0 | 2 | 2 | 3 | 4 |
| ・・・ | ・・・ | ・・・ | ・・・ | ・・・ | ・・・ | ・・・ | ・・・ |
| 36 | B12 | 疾患B | 0 | 7 | 7 | 5 | 4 |
Step 2. データを群に分ける
次に、読み込んだデータフレームを「正常」、「疾患A」、「疾患B」の 3 群に分けて、データフレーム配列 G_data に格納します。
# 読み込んだ検査データ(labo_data)のデータフレームを３群に分ける（G_j, j=1,2,3）
G_labels = ['正常', '疾患A', '疾患B']
G_data = []
for i, G_label in enumerate(G_labels):
G_data.append(T_data.query(f'群 == "{G_label}"').loc[:, '検査項目1':'検査項目5'])
Step 3. 射影軸（固有ベクトル）を求める
行列 $S_W^{-1} S_B$ の固有値と固有ベクトル（射影軸の方向を定めるベクトル）を求めます。
from pandas.core.window.rolling import window_aggregations
import numpy as np
# 配列の表示形式を設定：小数点以下2桁を出力する
np.set_printoptions(precision=2)
# グループ平均ベクトル G_mean[j]（\bar{x}^{(j)}）
G_mean = []
for i, data in enumerate(G_data):
G_mean.append(data.mean(axis='index'))
# 群サンプル数 n_data（\bar{n}）
n_data = np.array([len(df) for df in G_data])
# 全体平均ベクトル T_mean（\bar{x}）
T_mean = np.array([n_data[i] * mean.values for i, mean in enumerate(G_mean)]).sum(axis=0) / n_data.sum()
# 群間変動行列
Bs = []
for j in range(len(G_mean)):
v = G_mean[j].values -T_mean
B = n_data[j] * v.reshape(-1,1).dot(v.reshape(1,-1))
Bs.append(B)
B_matrix = sum(Bs)
# 群内変動行列
Ws = []
for j in range(len(G_mean)):
Ss = []
for i in range(len(G_data[j])):
v = G_data[j].values[i,:] - G_mean[j].values
S = v.reshape(-1,1).dot(v.reshape(1,-1))
Ss.append(S)
S_matrix = sum(Ss)
Ws.append(S_matrix)
W_matrix = sum(Ws)
# 行列（S_W^{-1}S_B）
W_matrix_inv = np.linalg.inv(W_matrix)
A = W_matrix_inv.dot(B_matrix)
# 行列（S_W^{-1}S_B）の固有値、固有ベクトルを求める関数
def eig_sort(A):
# A の固有値、固有ベクトルを求める（注意：np.linalg.eig は、固有ベクトルを列ベクトルとして返す）
evals, evecs = np.linalg.eig(A)
# 固有値の大きい順に並べ替える
idx = evals.argsort()[::-1]
# 固有ベクトルは行ベクトルで返す
return evals[idx], evecs.T[idx]
evals, evecs = eig_sort(A)
print(f'第1固有ベクトル：w1={evecs[0,:]}, 第1固有値：λ1={evals[0]:.3f}')
print(f'第2固有ベクトル：w1={evecs[1,:]}, 第2固有値：λ1={evals[1]:.3f}')
実行結果は次のようになります。
第1固有ベクトル：w1=[-0.39 -0.55 0.68 -0.23 0.14], 第1固有値：λ1=1.861
第2固有ベクトル：w1=[ 0.49 0.37 -0.06 -0.73 -0.29], 第2固有値：λ1=0.145
Step 4. サンプルをプロット
それでは、サンプルデータを、前節で求めた射影軸（固有ベクトル $\boldsymbol{w}_1, \boldsymbol{w}_2$）で構成される座標平面にプロットしてみましょう。
最初に、グラフ表示で日本語の文字化けを防ぐために、ライブラリ japanize-matplotlib をインストールします。
pip install japanize-matplotlib
Google Colaboratory を利用している場合は、上記のコマンドの先頭に "!" を付けると、コードセル内でインストールできます。
サンプルをプロットするコードは以下のとおりです。
import matplotlib.pyplot as plt
import japanize_matplotlib
from scipy.stats import multivariate_normal
# メッシュ・グリッド(xgrid,ygrid)に対応する2次元正規分布の確率密度を求める
def get_z(mean, cov):
meshgrid_z = []
for xrow, yrow in zip(meshgrid_x, meshgrid_y):
zrow = []
for xgrid, ygrid in zip(xrow,yrow):
zgrid = multivariate_normal(mean, cov).pdf(np.array([xgrid, ygrid]))
zrow.append(zgrid)
meshgrid_z.append(zrow)
meshgrid_z = np.array(meshgrid_z)
return meshgrid_z
# サンプルを２次元に射影し、２次元平面にサンプル (w_1^T x, w_2^T x) をプロット
# 固有ベクトル w_1, w_2
w = [evecs[0,:].reshape(1,-1), evecs[1,:].reshape(1,-1)]
# メッシュ・グリッド
w1s, w1e = -2.5, 2.5
w2s, w2e = -8, 2
meshgrid_x = np.linspace(w1s, w1e, 100)
meshgrid_y = np.linspace(w2s, w2e, 100)
meshgrid_x, meshgrid_y = np.meshgrid(meshgrid_x, meshgrid_y)
w1se = w1e - w1s
w2se = w2e - w2s
fsize_y = 6
fsize_x = fsize_y *w1se / w2se
plt.figure(figsize=(fsize_x, fsize_y))
# 群ごとに２次元平面にサンプルをプロット
for j in range(len(G_data)):
# サンプルの射影を散布図に描く
w1x = w[0].dot(G_data[j].values.T)
w2x = w[1].dot(G_data[j].values.T)
plt.scatter(w1x[0], w2x[0], label=G_labels[j])
# 等高線を描画
w12x = np.array([w1x[0], w2x[0]])
cov = np.cov(w12x)
mean = np.mean(w12x.T, axis=0)
meshgrid_z = get_z(mean, cov)
plt.contour(meshgrid_x, meshgrid_y, meshgrid_z, levels=5)
plt.xlabel(f'w_1^T x ({100 * evals[0] / evals.sum():.1f}%)')
plt.ylabel(f'w_2^T x ({100 * evals[1] / evals.sum():.1f}%)')
plt.title('正準判別分析')
plt.legend(loc="lower left")
結果を下図に示します。
図の横軸は固有ベクトル $\boldsymbol{w}_1$ 方向の射影軸で、縦軸は固有ベクトル $\boldsymbol{w}_2$ 方向の射影軸です。軸ラベルの丸がっこ内の数値は正準スコアの寄与率で、次式で定義されています。
$$
第 k 正準スコアの寄与率 = \frac{\lambda_k}{\sum_{k'=1}^{d} \lambda_{k'}}
\quad・・・\quad(22).
$$
正準スコアの寄与率は、各正準変数が元の変数にどれくらいの情報または変動を含んでいるかを示す指標です。図から、第 1 正準スコアの寄与率は $92.8\%$、第 2 正準スコアの寄与率は $7.2\%$で、両者を併せて $100\%$ となり、この２つの射影軸で全ての情報量を含んでいることが分かります。
図から、第 1 射影軸で「正常」群と「疾患A」群を判別し、第 2 射影軸で「正常」群あるいは「疾患A」群と「疾患B」群を判別できることが分かります。ただし、第 2 射影軸の正準スコアの寄与率は $7.2\%$ と小さいため、分解能が十分でないことがわかります。
23 判別分析（その１）へ 
参考文献
- Sugimoto Norio Art Production：18.2 重判別分析. http://www.snap-tck.com/room04/c01/stat/stat18/stat1802.html (2023.9.14 確認)
- 松井勇佑：やんわり理解する行列のランク. https://www.ite.or.jp/contents/kirari/2301kirari.pdf (2023.9.14 確認)
- おっぱいそん！：行列のランクと固有値の関係. https://oppython.hatenablog.com/entry/2015/10/23/175522 (2023.9.15 確認)
- Sugimoto Norio Art Production：第18章 重判別分析. http://www.snap-tck.com/room04/c01/stat/stat18/stat1801.html (2023.9.15 確認)