13. 主成分分析 (2)#

!pip install japanize-matplotlib
Requirement already satisfied: japanize-matplotlib in /usr/local/lib/python3.8/dist-packages (1.1.3)
Requirement already satisfied: matplotlib in /usr/local/lib/python3.8/dist-packages (from japanize-matplotlib) (3.4.2)
Requirement already satisfied: pyparsing>=2.2.1 in /usr/local/lib/python3.8/dist-packages (from matplotlib->japanize-matplotlib) (2.4.7)
Requirement already satisfied: cycler>=0.10 in /usr/local/lib/python3.8/dist-packages (from matplotlib->japanize-matplotlib) (0.10.0)
Requirement already satisfied: kiwisolver>=1.0.1 in /usr/local/lib/python3.8/dist-packages (from matplotlib->japanize-matplotlib) (1.3.1)
Requirement already satisfied: numpy>=1.16 in /home/okazaki/.local/lib/python3.8/site-packages (from matplotlib->japanize-matplotlib) (1.19.5)
Requirement already satisfied: pillow>=6.2.0 in /usr/local/lib/python3.8/dist-packages (from matplotlib->japanize-matplotlib) (8.3.1)
Requirement already satisfied: python-dateutil>=2.7 in /usr/local/lib/python3.8/dist-packages (from matplotlib->japanize-matplotlib) (2.8.1)
Requirement already satisfied: six in /home/okazaki/.local/lib/python3.8/site-packages (from cycler>=0.10->matplotlib->japanize-matplotlib) (1.15.0)
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import japanize_matplotlib

2次元のデータの例で主成分分析の基本的な考え方を掴んだところで、多次元のデータに対する主成分分析を考えよう。

13.1. データの表現#

\(d\)次元のベクトル\(\pmb{x} \in \mathbb{R}^d\)からなる事例を\(N\)件集めたデータを行列\(\pmb{X} \in \mathbb{R}^{N \times d}\)で表す。

(13.1)#\[\begin{align} \pmb{X} = \begin{pmatrix}\pmb{x}_1^\top \\ \pmb{x}_2^\top \\ \dots \\ \pmb{x}_N^\top \end{pmatrix} \end{align}\]

また、データの平均ベクトル\(\pmb{\bar{x}} \in \mathbb{R}^d\)を次式で定義する。

(13.2)#\[\begin{align} \pmb{\bar{x}} = \frac{1}{N} \sum_{i=1}^N \pmb{x}_i \end{align}\]

13.1.1. 射影したデータの分散#

あるベクトル\(\pmb{u} \in \mathbb{R}^d\)でデータ点を射影したとき、\(\pmb{u}\)上に射影された点の分散\(J\)は、

(13.3)#\[\begin{split} \begin{align} J &= \frac{1}{N} \sum_{i=1}^N \left(\pmb{u}^\top \pmb{x}_i - \pmb{u}^\top \pmb{\bar{x}}\right)^2 \\ &= \frac{1}{N} \sum_{i=1}^N \left\{\pmb{u}^\top (\pmb{x}_i - \pmb{\bar{x}})\right\}^2 \\ &= \frac{1}{N} \sum_{i=1}^N \left\{\pmb{u}^\top (\pmb{x}_i - \pmb{\bar{x}})\right\} \cdot \left\{\pmb{u}^\top (\pmb{x}_i - \pmb{\bar{x}})\right\} \\ &= \frac{1}{N} \sum_{i=1}^N \left\{\pmb{u}^\top (\pmb{x}_i - \pmb{\bar{x}})\right\} \cdot \left\{(\pmb{x}_i - \pmb{\bar{x}})^\top \pmb{u}\right\} \\ &= \pmb{u}^\top \left[\frac{1}{N} \sum_{i=1}^N (\pmb{x}_i - \pmb{\bar{x}})(\pmb{x}_i - \pmb{\bar{x}})^\top\right] \pmb{u} \\ &= \pmb{u}^\top \pmb{S} \pmb{u} \end{align} \end{split}\]

射影後の点の分散

\(N\)個のデータ点\(\pmb{x}_1, \dots, \pmb{x}_N \in \mathbb{R}^d\)をベクトル\(\pmb{u} \in \mathbb{R}^d\)で射影したとき、\(\pmb{u}\)上に射影された点の分散は、

(13.4)#\[ \begin{align} \pmb{u}^\top \pmb{S} \pmb{u} \end{align} \]

ここで、\(\pmb{S} \in \mathbb{R}^{d \times d}\)はデータの分散共分散行列である。

(13.5)#\[\begin{split} \begin{align} \pmb{S} &= \frac{1}{N} \sum_{i=1}^N (\pmb{x}_i - \pmb{\bar{x}}) (\pmb{x}_i - \pmb{\bar{x}})^\top \\ \end{align} \end{split}\]

なお、分散共分散行列\(\pmb{S}\)の各要素\(S_{a,b}\)は、行列\(\pmb{X}\)を用いると、

(13.6)#\[\begin{gather} S_{a,b} = \frac{1}{N} \sum_{i=1}^N (X_{i,a} - \bar{x}_a) (X_{i,b} - \bar{x}_b) \\ \bar{x}_a = \frac{1}{N} \sum_{i=1}^N X_{i,a}, \; \bar{x}_b = \frac{1}{N} \sum_{i=1}^N X_{i,b} \end{gather}\]

と表される。\(a\)\(b\)を入れ替えても同じ値が求まるので、分散共分散行列\(\pmb{S}\)\(d \times d\)の実対称行列である。

13.2. 固有値問題と主成分分析#

13.2.1. 第1主成分#

主成分分析では、ある単位ベクトル\(\pmb{u} \in \mathbb{R}^d\)でデータ点を射影したとき、射影された点の分散を最大化したい。前章で説明した通り、第1主成分\(\pmb{u}_1\)は、\(\|\pmb{u}_1\|=1\)の制約のもとで次式の目的関数を最大化する問題の解である。

(13.7)#\[\begin{align} J = \pmb{u}_1^\top \pmb{S} \pmb{u}_1 \end{align}\]

これを制約なしの最大化問題に書き換えるため、ラグランジュの未定乗数法を用いる。ラグランジュ関数\(\mathcal{L}\)は、

(13.8)#\[\begin{align} \mathcal{L} = \pmb{u}_1^\top \pmb{S} \pmb{u}_1 - \lambda_1 (\pmb{u}_1^\top\pmb{u}_1 - 1) \end{align}\]

ここで、\(\lambda_1\)はラグランジュ乗数である。ラグランジュ関数\(\mathcal{L}\)を最大化するため、\(\pmb{u}_1\)に関する偏微分を\(0\)とおくと、

(13.9)#\[\begin{align} \frac{\partial \mathcal{L}}{\partial \pmb{u}_1} = 2\pmb{S} \pmb{u}_1 - 2\lambda_1 \pmb{u}_1 &= 0 \\ \pmb{S} \pmb{u}_1 = \lambda_1 \pmb{u}_1 \end{align}\]

ゆえに、第1主成分を求める問題は分散共分散行列\(\pmb{S}\)の固有値問題に帰着する。なお、上式の両辺に左から\(\pmb{u}_1^\top\)をかけると、左辺は最大化したい目的関数\(J\)となり、右辺は\(\lambda_1\)となる。

(13.10)#\[ \begin{align} J = \pmb{u}_1^\top \pmb{S} \pmb{u}_1 = \pmb{u}_1^\top (\pmb{S} \pmb{u}_1) = \pmb{u}_1^\top (\lambda_1 \pmb{u}_1) = \lambda_1 \pmb{u}_1^\top \pmb{u}_1 = \lambda_1 \end{align} \]

したがって、射影後の分散\(J\)を最大にする第1主成分ベクトル\(\pmb{u}_1\)は、\(\pmb{S}\)の固有値の中で最大のものを\(\lambda_1\)として、それに対応する固有ベクトルを長さ\(1\)になるように正規化したものである。また、データ\(\pmb{X}\)を第1主成分ベクトル\(\pmb{u}_1\)に射影したとき、その軸における分散は\(\lambda_1\)である。

これまでの議論をまとめる。

主成分分析の第1主成分

\(N\)個の事例からなるデータ\(\pmb{X} = \begin{pmatrix}\pmb{x}_1 & \pmb{x}_2 & \dots & \pmb{x}_N \end{pmatrix}^\top\)が与えられ、その分散共分散行列を\(\pmb{S}\)とする。主成分分析の第1主成分ベクトル\(\pmb{u}_1\)は以下の固有値問題の解である。

(13.11)#\[ \begin{align} \pmb{S} \pmb{u}_1 = \lambda_1 \pmb{u}_1 \end{align} \]

ここで、\(\lambda_1\)\(\pmb{S}\)の最大の固有値で、\(\pmb{u}_1\)\(\lambda_1\)に対応する固有ベクトルを正規化したものである。データの平均ベクトルを\(\bar{\pmb{x}}\)と書くことにすると、データ\(\pmb{X}\)に含まれるある事例\(\pmb{x}\)の第1主成分得点は、

(13.12)#\[ \begin{align} \pmb{u}_1^\top (\pmb{x} - \pmb{\bar{x}}) \end{align} \]

データ点を第1主成分ベクトル\(\pmb{u}_1\)に射影したとき、その軸におけるデータの分散は\(\lambda_1\)である。

13.2.2. 第2主成分以降#

第2主成分以降を一般的に求めるため、第\(k\)主成分ベクトル\(\pmb{u}_1, \dots, \pmb{u}_k\)までを求めてあると仮定し、第\(k+1\)主成分ベクトルを求める(\(1 \leq k < d\))。各主成分ベクトルの長さは\(1\)で、互いに直交するように選ぶので、主成分ベクトルは以下の制約を満たす(つまり、\(\pmb{u}_1, \dots, \pmb{u}_k\)は正規直交系である)。

(13.13)#\[\begin{split} \begin{align} \pmb{u}^\top_i \pmb{u}_j = \begin{cases} 1 & (i = j) \\ 0 & (i \neq j) \end{cases} \end{align} \end{split}\]

\(k+1\)主成分ベクトル\(\pmb{u}_{k+1}\)は、\(\|\pmb{u}_{k+1}\|=1, \pmb{u}_1^\top \pmb{u}_{k+1} = \dots = \pmb{u}_k^\top \pmb{u}_{k+1} = 0\)の制約のもとで次の目的関数を最大化する問題の解である。

(13.14)#\[\begin{align} J = \pmb{u}_{k+1}^\top \pmb{S} \pmb{u}_{k+1} \end{align}\]

これを制約なしの最大化問題に書き換えるため、ラグランジュの未定乗数法を用いる。ラグランジュ関数\(\mathcal{L}\)は、

(13.15)#\[\begin{align} \mathcal{L} = \pmb{u}_{k+1}^\top \pmb{S} \pmb{u}_{k+1} - \lambda_{k+1} (\pmb{u}_{k+1}^\top\pmb{u}_{k+1} - 1) - \sum_{i=1}^k \alpha_i \pmb{u}_i^\top \pmb{u}_{k+1} \end{align}\]

ただし、\(\alpha_1, \dots, \alpha_k, \lambda_{k+1}\)はラグランジュ乗数である。ラグランジュ関数\(\mathcal{L}\)を最大化するため、\(\pmb{u}_{k+1}\)に関する偏微分を\(0\)とおくと、

(13.16)#\[ \begin{align} \frac{\partial \mathcal{L}}{\partial \pmb{u}_{k+1}} = 2\pmb{S} \pmb{u}_{k+1} - 2\lambda_{k+1} \pmb{u}_{k+1} - \sum_{i=1}^k \alpha_i \pmb{u}_i &= 0 \end{align} \]

ここで、上の等式において\(\pmb{u}_1^\top\)を左からかけ、式(13.13)に注意しながら展開すると、

(13.17)#\[\begin{align} 2\pmb{u}_1^\top \pmb{S} \pmb{u}_{k+1} - 2\lambda_{k+1} \pmb{u}_1^\top \pmb{u}_{k+1} - \sum_{i=1}^k \alpha_i \pmb{u}_1^\top \pmb{u}_i &= 0 \\ 2\pmb{u}_1^\top \pmb{S} \pmb{u}_{k+1} - \alpha_1 &= 0 \\ \alpha_1 &= 2\pmb{u}_1^\top \pmb{S} \pmb{u}_{k+1} = 2\pmb{u}_{k+1}^\top \pmb{S} \pmb{u}_1 = 2\pmb{u}_{k+1}^\top \lambda_1 \pmb{u}_1 = 0 \end{align}\]

同様に、全ての\(j \in \{1, \dots, k\}\)に対して、\(\pmb{u}_j^\top\)を左からかけることで、\(\alpha_j = 0\)が示される。ゆえに、式(13.16)は、

(13.18)#\[\begin{align} \pmb{S} \pmb{u}_{k+1} = \lambda_{k+1} \pmb{u}_{k+1} \end{align}\]

と整理される。したがって、第\(k+1\)主成分以降を求める問題も分散共分散行列\(\pmb{S}\)の固有値問題に帰着する。これを\(k=1\)から\(k=d-1\)まで繰り返すことで、すべての主成分ベクトルを求めることができる。

固有値問題による主成分分析をまとめる。

主成分分析と固有値問題

\(N\)個の事例からなるデータ\(\pmb{X} = \begin{pmatrix}\pmb{x}_1 & \pmb{x}_2 & \dots & \pmb{x}_N \end{pmatrix}^\top\)が与えられ、その分散共分散行列を\(\pmb{S}\)とする。\(\pmb{S}\)\(d\)個の固有値を大きい順に並べ、

(13.19)#\[ \begin{align} \lambda_1 \geq \lambda_2 \geq \dots \geq \lambda_d \geq 0 \end{align} \]

とする。また、\(\lambda_1, \lambda_2, \dots, \lambda_d\)に対応する固有ベクトルを\(\pmb{u}_1, \pmb{u}_1, \dots, \pmb{u}_d\)とする。ただし、固有ベクトルは、

(13.20)#\[\begin{align} \pmb{u}^\top_i \pmb{u}_j = \begin{cases} 1 & (i = j) \\ 0 & (i \neq j) \end{cases} \end{align}\]

を満たすように求める。すると、\(\pmb{u}_1, \pmb{u}_1, \dots, \pmb{u}_d\)はそれぞれ、第1主成分、第2主成分、……、第\(d\)主成分に対応する。データの平均ベクトルを\(\bar{\pmb{x}}\)と書くことにすると、データ\(\pmb{X}\)に含まれるある事例\(\pmb{x}\)の第\(i\)主成分得点は(\(i \in \{1, \dots, d\}\))、

(13.21)#\[ \begin{align} \pmb{u}_i^\top (\pmb{x} - \pmb{\bar{x}}) \end{align} \]

データ\(\pmb{X}\)を各主成分に射影したとき、各軸におけるデータの分散はそれぞれ、\(\lambda_1, \lambda_2, \dots, \lambda_d\)である。

なお、式(13.19)において、\(\pmb{S}\)の固有値が実数で、かつ全て非負であることは後ほど説明する。

13.3. スペクトル分解と主成分分析#

主成分分析への理解を深めるため、スペクトル分解との関連を説明する。

13.3.1. スペクトル分解#

まず、スペクトル分解を復習する。一般に、実対象行列は以下の性質を持つ。

ある実対称行列\(\pmb{S} \in \mathbb{R}^{d \times d}\)の固有値を\(\lambda_1, \lambda_2, \dots, \lambda_d\)とし、対応する固有ベクトルを\(\pmb{q}_1, \pmb{q}_2, \dots, \pmb{q}_d \in \mathbb{R}^d\)とする。

(13.22)#\[ \begin{align} \pmb{S} \pmb{q}_1 = \lambda_1 \pmb{q}_1, \; \pmb{S} \pmb{q}_2 = \lambda_2 \pmb{q}_2, \; \dots, \; \pmb{S} \pmb{q}_d = \lambda_d \pmb{q}_d \end{align} \]

なお、固有ベクトルを正規直交系を構成するように選ぶことができるので、\(i, j \in \{1, \dots, d\}\)に関して、

(13.23)#\[\begin{split} \begin{align} \pmb{q}_i^\top\pmb{q}_j = \begin{cases} 1 & (i = j) \\ 0 & (i \neq j) \end{cases} \end{align} \end{split}\]

等式(13.22)をベクトルとして横に並べると、次のように整理できる。

(13.24)#\[\begin{split} \begin{align} \begin{pmatrix} \pmb{S} \pmb{q}_1 & \pmb{S} \pmb{q}_2 & \dots & \pmb{S} \pmb{q}_d \end{pmatrix} &= \begin{pmatrix} \lambda_1 \pmb{q}_1 & \lambda_2 \pmb{q}_2 & \dots & \lambda_d \pmb{q}_d \end{pmatrix} \\ \pmb{S} \begin{pmatrix} \pmb{q}_1 & \pmb{q}_2 & \dots & \pmb{q}_d \end{pmatrix} &= \begin{pmatrix} \pmb{q}_1 & \pmb{q}_2 & \dots & \pmb{q}_d \end{pmatrix} \begin{pmatrix} \lambda_1 & 0 & \dots & 0 \\ 0 & \lambda_2 & \dots & 0 \\ 0 & 0 & \dots & 0 \\ 0 & 0 & \dots & \lambda_d \\ \end{pmatrix} \\ \pmb{S} \pmb{Q} &= \pmb{Q} \pmb{\Lambda} \end{align} \end{split}\]

ここで、\(\pmb{\Lambda} \in \mathbb{R}^{d \times d}\)は固有値を要素とする対角行列である。

(13.25)#\[\begin{align} \pmb{\Lambda} = {\rm diag}(\lambda_1, \lambda_2, \dots, \lambda_d) = \begin{pmatrix} \lambda_1 & 0 & \ldots & 0 \\ 0 & \lambda_2 & \ldots & 0 \\ \vdots & \vdots & \ddots & \vdots \\ 0 & 0 & \ldots & \lambda_d \\ \end{pmatrix} \end{align}\]

また、行列\(\pmb{Q} \in \mathbb{R}^{d \times d}\)は列ベクトルである固有ベクトルを横方向に配置したものである。

(13.26)#\[\begin{align} \pmb{Q} = \begin{pmatrix} \pmb{q}_1 & \pmb{q}_2 & \dots & \pmb{q}_d \end{pmatrix} \end{align}\]

以下に示すように、\(\pmb{Q}^\top \pmb{Q} = \pmb{I}\)であるから、\(\pmb{Q}\)は直交行列である。

(13.27)#\[\begin{split} \begin{align} \pmb{Q}^\top \pmb{Q} = \begin{pmatrix} \pmb{q}_1^\top \\ \pmb{q}_2^\top \\ \dots \\ \pmb{q}_d^\top \end{pmatrix} \begin{pmatrix} \pmb{q}_1 & \pmb{q}_2 & \dots & \pmb{q}_d \end{pmatrix} = \begin{pmatrix} \pmb{q}_1^\top\pmb{q}_1 & \pmb{q}_1^\top\pmb{q}_2 & \ldots & \pmb{q}_1^\top\pmb{q}_d \\ \pmb{q}_2^\top\pmb{q}_1 & \pmb{q}_2^\top\pmb{q}_2 & \ldots & \pmb{q}_2^\top\pmb{q}_d \\ \vdots & \vdots & \ddots & \vdots \\ \pmb{q}_d^\top\pmb{q}_1 & \pmb{q}_d^\top\pmb{q}_2 & \ldots & \pmb{q}_d^\top\pmb{q}_d \\ \end{pmatrix} = \begin{pmatrix} 1 & 0 & \ldots & 0 \\ 0 & 1 & \ldots & 0 \\ \vdots & \vdots & \ddots & \vdots \\ 0 & 0 & \ldots & 1 \\ \end{pmatrix} = \pmb{I} \end{align} \end{split}\]

(13.24)の右から\(\pmb{Q}^\top\)をかけ、左辺において\(\pmb{Q}\pmb{Q}^\top = \pmb{I}\)であることに注意すると、スペクトル分解(spectral decomposition)が得られる。

スペクトル分解

実対称行列\(\pmb{S} \in \mathbb{R}^{d \times d}\)は、直交行列\(\pmb{Q} \in \mathbb{R}^{d \times d}\)、対角行列\(\pmb{\Lambda} \in \mathbb{R}^{d \times d}\)、直交行列\(\pmb{Q}^\top\)の積に分解できる。

(13.28)#\[\begin{align} \pmb{S} &= \pmb{Q} \pmb{\Lambda} \pmb{Q}^\top \\ &= \begin{pmatrix} \vert & \vert & \ldots & \vert \\ \pmb{q}_1 & \pmb{q}_2 & \ldots & \pmb{q}_d \\ \vert & \vert & \ldots & \vert \\ \end{pmatrix} \begin{pmatrix} \lambda_1 & 0 & \ldots & 0 \\ 0 & \lambda_2 & \ldots & 0 \\ \vdots & \vdots & \ddots & \vdots \\ 0 & 0 & \ldots & \lambda_d \\ \end{pmatrix} \begin{pmatrix} \text{―} & \pmb{q}_1^\top & \text{―} \\ \text{―} & \pmb{q}_2^\top & \text{―} \\ \vdots & \vdots & \vdots \\ \text{―} & \pmb{q}_d^\top & \text{―} \end{pmatrix} \\ &= \lambda_1 \begin{pmatrix} \vert \\ \pmb{q}_1 \\ \vert \\ \end{pmatrix} \begin{pmatrix} \text{―} \; \pmb{q}_1^\top \; \text{―} \end{pmatrix} + \lambda_2 \begin{pmatrix} \vert \\ \pmb{q}_2 \\ \vert \\ \end{pmatrix} \begin{pmatrix} \text{―} \; \pmb{q}_2^\top \; \text{―} \end{pmatrix} + \dots + \lambda_d \begin{pmatrix} \vert \\ \pmb{q}_d \\ \vert \\ \end{pmatrix} \begin{pmatrix} \text{―} \; \pmb{q}_d^\top \; \text{―} \end{pmatrix} \\ &= \sum_{i = 1}^{d} \lambda_i \pmb{q}_i \pmb{q}_i^\top \end{align}\]

なお、

  • \(\pmb{\Lambda}\)の対角要素\(\lambda_1, \lambda_2, \dots, \lambda_d\)\(\pmb{S}\)の固有値を並べたものである。

  • \(\pmb{Q}\)の列ベクトル\(\pmb{q}_1, \dots, \pmb{q}_d \in \mathbb{R}^d\)\(\lambda_1, \lambda_2, \dots, \lambda_d\)に対応する固有ベクトルであり、正規直交系をなす。

13.3.2. 主成分分析の導出#

スペクトル分解に基づいて主成分分析を導出してみよう。主成分分析では、データ行列\(\pmb{X} \in \mathbb{R}^{N \times d}\)に対して式(13.5)で計算される分散共分散行列\(\pmb{S} \in \mathbb{R}^{d \times d}\)が出発点である。以下の二つの性質により、分散共分散行列\(\pmb{S}\)の固有値は非負である。

分散共分散行列\(\pmb{S}\)は実対称行列であるから、直交行列\(\pmb{Q} \in \mathbb{R}^{d \times d}\)と対角行列\(\pmb{\Lambda} = {\rm diag}(\lambda_1, \dots, \lambda_d) \in \mathbb{R}^{d \times d}\)によるスペクトル分解が存在する。

(13.32)#\[\begin{align} \pmb{S} = \pmb{Q} \pmb{\Lambda} \pmb{Q}^\top \end{align}\]

さらに、分散共分散行列\(\pmb{S}\)は半正定値対称行列であり、その固有値はすべて非負である。そこで、\(\pmb{S}\)の固有値を大きい順に並べておく。

(13.33)#\[ \begin{align} \lambda_1 \geq \lambda_2 \geq \dots \geq \lambda_d \geq 0 \end{align} \]

また、これらの固有値に対応するように固有ベクトル\(\pmb{Q} = \begin{pmatrix}\pmb{q}_1 & \pmb{q}_2 & \dots & \pmb{q}_d\end{pmatrix}\)を並び替えておく。

さて、式(13.38)で示したように、主成分分析における第1主成分\(\pmb{u}_1 \in \mathbb{R}^d\)は以下の最大化問題で求まる。

(13.34)#\[\begin{align} \pmb{u}_1 = \mathop{\rm argmax}\limits_{\|\pmb{u}\| = 1} \pmb{u}^\top \pmb{S} \pmb{u} = \mathop{\rm argmax}\limits_{\|\pmb{u}\| \neq 0} \frac{\pmb{u}^\top \pmb{S} \pmb{u}}{\pmb{u}^\top \pmb{u}} \end{align}\]

ここで、\({\rm argmax}\)の中身、

(13.35)#\[\begin{align} \frac{\pmb{u}^\top \pmb{S} \pmb{u}}{\pmb{u}^\top \pmb{u}} \end{align}\]

はレイリー商(Rayleigh quotient)と呼ばる。以下で示すように、このレイリー商は\(\pmb{u} = \pm \pmb{q}_1\)のとき最大値\(\lambda_1\)をとる。すなわち、データ\(\pmb{X}\)を第1主成分ベクトルに射影するとき、その分散の最大値は\(\lambda_1\)であり、その最大値を与える第1主成分ベクトルは\(\pmb{u}_1 = \pm \pmb{q}_1\)である。

レイリー商による最大値

(13.36)#\[\begin{align} \frac{\pmb{u}^\top \pmb{S} \pmb{u}}{\pmb{u}^\top \pmb{u}} &= \frac{\pmb{u}^\top (\pmb{Q}\pmb{\Lambda}\pmb{Q}^\top) \pmb{u}}{\pmb{u}^\top \pmb{u}} && (\because \pmb{S} = \pmb{Q}\pmb{\Lambda}\pmb{Q}^\top) \\ &= \frac{\pmb{u}^\top \pmb{Q}\pmb{\Lambda}\pmb{Q}^\top \pmb{u}}{\pmb{u}^\top\pmb{Q} \pmb{Q}^\top\pmb{u}} && (\because \pmb{Q} \pmb{Q}^\top = \pmb{I}) \\ &= \frac{(\pmb{Q}^\top \pmb{u})^\top\pmb{\Lambda}\pmb{Q}^\top \pmb{u}}{(\pmb{Q}^\top \pmb{u})^\top \pmb{Q}^\top\pmb{u}} && (\because \pmb{Q} \pmb{Q}^\top = \pmb{I}) \\ &= \frac{\pmb{z}^\top \pmb{\Lambda}\pmb{z}}{\pmb{z}^\top\pmb{z}} && (\pmb{z} := \pmb{Q}^\top \pmb{u}) \\ &= \frac{\lambda_1 z_1^2 + \lambda_2 z_2^2 + \dots + \lambda_d z_d^2}{z_1^2 + z_2^2 + \dots + z_d^2} && (\because \pmb{\Lambda} = {\rm diag}(\lambda_1, \lambda_2, \dots, \lambda_d)) \end{align}\]

上式は\(\lambda_1, \lambda_2, \dots, \lambda_d\)\(z_1^2, z_2^2, \dots, z_d^2\)による重み付き平均である。したがって、その最大値は、

(13.37)#\[\begin{align} \max_{\|\pmb{u}\| \neq 0} \frac{\pmb{u}^\top \pmb{S} \pmb{u}}{\pmb{u}^\top \pmb{u}} = \max_{\|\pmb{z}\| \neq 0} \frac{\lambda_1 z_1^2 + \lambda_2 z_2^2 + \dots + \lambda_d z_d^2}{z_1^2 + z_2^2 + \dots + z_d^2} \leq \lambda_1 \\ \end{align}\]

等号が成立するのは、\(z_1^2 = 1, z_2^2 = z_3^2 = \dots = z_d^2 = 0\)、すなわち\(\pmb{u} = \pm \pmb{q}_1\)のときである。

続いて、主成分分析で第\(k\)主成分までデータ\(\pmb{X}\)を射影する場合を考える(\(1 \leq k \leq d\))。\(k\)個の直交基底ベクトルで射影した後のデータの分散の和の最大値は\(\lambda_1 + \dots + \lambda_k\)であること、その最大値は\(\pmb{q}_1, \dots, \pmb{q}_k\)で射影した場合であることを示す [Dasgupta, 2008]

\(k\)個の直交基底ベクトルを\(\pmb{u}_1, \dots, \pmb{u}_k\)とおく。式(13.4)より、これらのベクトルでデータを射影した後の分散の和\(J_k\)は、分散共分散行列\(\pmb{S}\)のスペクトル分解\(\pmb{Q}\pmb{\Lambda}\pmb{Q}^\top\)を用いて、次のように展開できる。

(13.38)#\[\begin{split} \begin{align} J_k &= \sum_{i=1}^k \pmb{u}_i^\top \pmb{S} \pmb{u}_i \\ &= \sum_{i=1}^k \pmb{u}_i^\top (\pmb{Q} \pmb{\Lambda} \pmb{Q}^\top) \pmb{u}_i && (\because \pmb{S} = \pmb{Q}\pmb{\Lambda}\pmb{Q}^\top) \\ &= \sum_{i=1}^k \pmb{u}_i^\top \left( \sum_{j = 1}^{d} \lambda_j \pmb{q}_j \pmb{q}_j^\top \right) \pmb{u}_i && (\mbox{スペクトル分解の別の形}) \\ &= \sum_{j = 1}^{d} \lambda_j \sum_{i=1}^k (\pmb{u}_i^\top\pmb{q}_j)^2 && (\because \pmb{u}_i^\top\pmb{q}_j = \pmb{q}_j^\top\pmb{u}_i) \\ &= \sum_{j = 1}^{d} \lambda_j z_j && (z_j := \sum_{i=1}^k (\pmb{u}_i^\top\pmb{q}_j)^2) \\ \end{align} \end{split}\]

ゆえに、分散の和\(J_k\)\(d\)個の固有値\(\lambda_1, \dots, \lambda_d\)\(z_1, \dots, z_d\)による重み付き和として表される。そこで、\(j \in \{1, 2, \dots, d\}\)に関して\(z_j\)の性質を調べる。まず、ベッセルの不等式より、

(13.39)#\[\begin{align} z_j = \sum_{i=1}^k (\pmb{u}_i^\top\pmb{q}_j)^2 \leq \|\pmb{q}_j\|^2 = 1 \end{align}\]

さらに、\(z_j\)の和は、

(13.40)#\[\begin{align} \sum_{j = 1}^{d} z_j = \sum_{j = 1}^{d} \sum_{i=1}^k (\pmb{u}_i^\top\pmb{q}_j)^2 = \sum_{i=1}^k \sum_{j = 1}^{d} \pmb{u}_i^\top\pmb{q}_j \pmb{q}_j^\top\pmb{u}_i = \sum_{i=1}^k \pmb{u}_i^\top\pmb{Q} \pmb{Q}^\top\pmb{u}_i = \sum_{i=1}^k \pmb{u}_i^\top\pmb{u}_i = k \end{align}\]

したがって、式(13.38)\(0 \leq z_j \leq 1\)かつ、\(\sum_{j = 1}^{d} z_j = k\)による固有値\(\lambda_1, \dots, \lambda_d\)の重み付き和である。式(13.33)のように、分散共分散行列\(\pmb{S}\)の固有値を大きい順に並べてある場合、式(13.38)が最大となるのは、\(z_1 = \dots = z_k = 1, z_{k+1} = \dots = z_d = 0\)と選んだときである。\(1 \leq j \leq k\)に対して、\(z_j = 1\)を満たすには、\(\pmb{u}_j = \pmb{q}_j\)とすればよい。

まとめると、\(k\)個の直交基底ベクトルでデータ\(\pmb{X}\)を射影した後の分散の和\(J_k\)の最大値は、

(13.41)#\[ \begin{align} \max J_k = \sum_{i=1}^k \lambda_i \end{align} \]

この最大値は、データを\(\pmb{q}_1, \dots, \pmb{q}_k\)で射影した時に得られる。

固有値問題と主成分分析では、データの分散共分散行列\(\pmb{S}\)の固有値を大きい順に並べ、対応する固有ベクトルを主成分ベクトルとして一つずつ選んでいた。式(13.41)は、\(k\)個の直交基底ベクトルでデータを射影するとき、\(\pmb{q}_1, \pmb{q}_2, \dots, \pmb{q}_k\)を射影ベクトルに選ぶことで、射影後のデータの分散の和が最大となること、その分散の和の最大値が\(\lambda_1 + \lambda_2 + \dots + \lambda_k\)であることを保証している。

13.4. 主成分分析と次元圧縮#

データ\(\pmb{X}\)の共分散行列\(\pmb{S}\)の対角成分の和は、データ\(\pmb{X}\)の分散の和である。\(\pmb{S} = \pmb{Q} \pmb{\Lambda} \pmb{Q}^\top\)のスペクトル分解を用いると、\(\pmb{S}\)の対角成分の和は、

(13.42)#\[\begin{align} \mathrm{tr}(\pmb{S}) &= \mathrm{tr}(\pmb{Q} \pmb{\Lambda} \pmb{Q}^\top) \\ &= \mathrm{tr}(\lambda_1 \pmb{q}_1 \pmb{q}_1^\top + \lambda_2 \pmb{q}_2 \pmb{q}_2^\top + \dots + \lambda_d \pmb{q}_d \pmb{q}_d^\top) \\ &= \lambda_1 \pmb{q}_1^\top \pmb{q}_1 + \lambda_2 \pmb{q}_2^\top \pmb{q}_2 + \dots + \lambda_d \pmb{q}_d^\top \pmb{q}_d \\ &= \lambda_1 + \lambda_2 + \dots + \lambda_d \\ \end{align}\]

つまり、データ\(\pmb{X}\)の分散の和は分散共分散行列\(\pmb{S}\)の固有値の和に等しい

\(k\)主成分に対応する固有値\(\lambda_k\)は、元のデータ\(X\)を第\(k\)主成分の軸に射影したときの分散を表す。全ての主成分における分散の総和に対し、第\(k\)主成分の分散が占める割合を第\(k\)主成分の寄与率と呼ぶ。

寄与率

(13.43)#\[\begin{align} \frac{\lambda_k}{\sum_{i=1}^d \lambda_i} = \frac{\lambda_k}{\lambda_1 + \lambda_2 + \dots + \lambda_d} \end{align}\]

また、第\(k\)主成分までに表現された分散の総和が、データ全体の分散の総和に占める割合を第\(k\)主成分までの累積寄与率と呼ぶ。

累積寄与率

(13.44)#\[\begin{align} \frac{\sum_{i=1}^k \lambda_i}{\sum_{i=1}^d \lambda_i} = \frac{\lambda_1 + \lambda_2 \dots + \lambda_k}{\lambda_1 + \lambda_2 + \dots + \lambda_d} \end{align}\]

累積寄与率は、元々\(d\)次元のベクトルで表されていたデータ点を主成分分析で\(k\)次元に圧縮したとき、元データの分散をどの程度表現できているかを表す指標である。

13.5. 応用例 (1): 麺類の支出額#

前章の冒頭で紹介した教育用標準データセット(SSDSE: Standardized Statistical Data Set for Education)の「C. 都道府県庁所在市別、家計消費データ」から麺類だけを抜き出したデータを2次元平面上で可視化してみよう。まず、データを読み込んで表示する。

df = pd.read_excel('https://www.nstac.go.jp/SSDSE/data/2021/SSDSE-C-2021.xlsx', skiprows=1)
df = df[['都道府県', '生うどん・そば', "乾うどん・そば", "パスタ", "中華麺", "カップ麺", "即席麺", "他の麺類"]]
df = df.iloc[1:]
df
都道府県 生うどん・そば 乾うどん・そば パスタ 中華麺 カップ麺 即席麺 他の麺類
1 北海道 3162 2082 1266 4152 5189 1609 726
2 青森県 2964 2224 1114 5271 6088 2039 554
3 岩手県 3349 2475 1305 5991 5985 1889 898
4 宮城県 3068 2407 1339 5017 5295 1709 842
5 秋田県 3231 3409 1019 5172 5013 1680 554
6 山形県 4478 3084 1288 5236 5875 1745 754
7 福島県 2963 2705 1064 4397 5862 1687 919
8 茨城県 3353 2477 1248 4034 4562 1440 653
9 栃木県 3908 2218 1391 4534 4945 1860 742
10 群馬県 4563 1948 1203 4153 5049 1544 546
11 埼玉県 4016 2256 1487 4512 4584 1568 949
12 千葉県 3389 2277 1441 4582 4513 1840 835
13 東京都 3088 2385 1595 4592 3953 1734 974
14 神奈川県 3410 2684 1472 4484 4035 1741 1049
15 新潟県 2697 3409 1308 4192 6095 2201 612
16 富山県 3294 2738 1184 4367 5576 1928 491
17 石川県 3597 2233 1339 5159 4408 1880 607
18 福井県 4124 2549 1200 4247 3838 1875 387
19 山梨県 3870 1947 1123 4529 4786 1580 660
20 長野県 4075 2705 1267 4333 4346 1620 808
21 岐阜県 3154 2106 1239 3816 4566 1907 650
22 静岡県 3045 2414 1319 5106 4582 1682 935
23 愛知県 3991 2198 1394 4239 4490 1882 731
24 三重県 3720 2190 1251 4176 4222 1922 595
25 滋賀県 3835 2004 1295 4121 4397 1976 659
26 京都府 3816 2160 1345 4136 3216 1959 596
27 大阪府 3999 1599 1013 3567 5564 2246 540
28 兵庫県 3437 1988 1129 3914 4630 1984 795
29 奈良県 3543 2727 1229 3973 4067 2003 617
30 和歌山県 3024 2363 940 3603 3972 1921 502
31 鳥取県 3455 1684 1251 3700 5687 3007 436
32 島根県 3456 2314 1316 4166 4233 2086 421
33 岡山県 3384 2219 1425 3958 4655 2036 522
34 広島県 3285 1851 1404 4036 4047 2077 585
35 山口県 3426 1840 1164 3342 4937 2533 535
36 徳島県 3692 2106 1107 3340 4149 2038 409
37 香川県 6558 3940 972 3661 4412 2170 411
38 愛媛県 3744 1939 1197 3624 3728 1825 498
39 高知県 3333 1513 1134 3306 5126 2234 478
40 福岡県 2877 1745 1354 3376 4084 2128 640
41 佐賀県 2428 2032 1132 3251 4329 2330 517
42 長崎県 2556 1882 1066 3720 4054 1893 463
43 熊本県 2307 1647 1303 2840 4631 2864 501
44 大分県 2493 1597 1074 3164 4158 2384 552
45 宮崎県 2668 1573 1185 3290 4098 2232 463
46 鹿児島県 2868 2391 1289 3524 3994 2054 462
47 沖縄県 1753 1488 1012 4424 3342 1481 229

7種類の麺類の品目に関して、47都道府県庁所在市別の年間支出金額が収録されたデータである。したがって、データ\(\pmb{X}\)\(47 \times 7\)の行列で表される。

X = df.iloc[:,1:].values
X
array([[3162, 2082, 1266, 4152, 5189, 1609,  726],
       [2964, 2224, 1114, 5271, 6088, 2039,  554],
       [3349, 2475, 1305, 5991, 5985, 1889,  898],
       [3068, 2407, 1339, 5017, 5295, 1709,  842],
       [3231, 3409, 1019, 5172, 5013, 1680,  554],
       [4478, 3084, 1288, 5236, 5875, 1745,  754],
       [2963, 2705, 1064, 4397, 5862, 1687,  919],
       [3353, 2477, 1248, 4034, 4562, 1440,  653],
       [3908, 2218, 1391, 4534, 4945, 1860,  742],
       [4563, 1948, 1203, 4153, 5049, 1544,  546],
       [4016, 2256, 1487, 4512, 4584, 1568,  949],
       [3389, 2277, 1441, 4582, 4513, 1840,  835],
       [3088, 2385, 1595, 4592, 3953, 1734,  974],
       [3410, 2684, 1472, 4484, 4035, 1741, 1049],
       [2697, 3409, 1308, 4192, 6095, 2201,  612],
       [3294, 2738, 1184, 4367, 5576, 1928,  491],
       [3597, 2233, 1339, 5159, 4408, 1880,  607],
       [4124, 2549, 1200, 4247, 3838, 1875,  387],
       [3870, 1947, 1123, 4529, 4786, 1580,  660],
       [4075, 2705, 1267, 4333, 4346, 1620,  808],
       [3154, 2106, 1239, 3816, 4566, 1907,  650],
       [3045, 2414, 1319, 5106, 4582, 1682,  935],
       [3991, 2198, 1394, 4239, 4490, 1882,  731],
       [3720, 2190, 1251, 4176, 4222, 1922,  595],
       [3835, 2004, 1295, 4121, 4397, 1976,  659],
       [3816, 2160, 1345, 4136, 3216, 1959,  596],
       [3999, 1599, 1013, 3567, 5564, 2246,  540],
       [3437, 1988, 1129, 3914, 4630, 1984,  795],
       [3543, 2727, 1229, 3973, 4067, 2003,  617],
       [3024, 2363,  940, 3603, 3972, 1921,  502],
       [3455, 1684, 1251, 3700, 5687, 3007,  436],
       [3456, 2314, 1316, 4166, 4233, 2086,  421],
       [3384, 2219, 1425, 3958, 4655, 2036,  522],
       [3285, 1851, 1404, 4036, 4047, 2077,  585],
       [3426, 1840, 1164, 3342, 4937, 2533,  535],
       [3692, 2106, 1107, 3340, 4149, 2038,  409],
       [6558, 3940,  972, 3661, 4412, 2170,  411],
       [3744, 1939, 1197, 3624, 3728, 1825,  498],
       [3333, 1513, 1134, 3306, 5126, 2234,  478],
       [2877, 1745, 1354, 3376, 4084, 2128,  640],
       [2428, 2032, 1132, 3251, 4329, 2330,  517],
       [2556, 1882, 1066, 3720, 4054, 1893,  463],
       [2307, 1647, 1303, 2840, 4631, 2864,  501],
       [2493, 1597, 1074, 3164, 4158, 2384,  552],
       [2668, 1573, 1185, 3290, 4098, 2232,  463],
       [2868, 2391, 1289, 3524, 3994, 2054,  462],
       [1753, 1488, 1012, 4424, 3342, 1481,  229]])

\(\pmb{X}\)の各事例に対応する都道府県名を変数Cに保存しておく。

C = df.iloc[:,0].values
C
array(['北海道', '青森県', '岩手県', '宮城県', '秋田県', '山形県', '福島県', '茨城県', '栃木県',
       '群馬県', '埼玉県', '千葉県', '東京都', '神奈川県', '新潟県', '富山県', '石川県', '福井県',
       '山梨県', '長野県', '岐阜県', '静岡県', '愛知県', '三重県', '滋賀県', '京都府', '大阪府',
       '兵庫県', '奈良県', '和歌山県', '鳥取県', '島根県', '岡山県', '広島県', '山口県', '徳島県',
       '香川県', '愛媛県', '高知県', '福岡県', '佐賀県', '長崎県', '熊本県', '大分県', '宮崎県',
       '鹿児島県', '沖縄県'], dtype=object)

\(\pmb{X}\)の品目名を変数Fに保存しておく。

F = df.columns[1:].values
F
array(['生うどん・そば', '乾うどん・そば', 'パスタ', '中華麺', 'カップ麺', '即席麺', '他の麺類'],
      dtype=object)

データの傾向を掴むために、散布図行列を描画してみる。対角部分はヒストグラム、非対角部分は2つの変数間に着目した場合の散布図が示される。

from pandas import plotting 
plotting.scatter_matrix(df.iloc[:,1:], figsize=(10, 10), alpha=0.5)
plt.show()
../_images/04pca2_21_0.png

sklearn.decomposition.PCAを用いて主成分分析を実行する。

from sklearn.decomposition import PCA

pca = PCA()
P = pca.fit_transform(X)

7個の主成分ベクトル(PC1からPC7)を表示する。第1主成分ベクトルの係数が概ね正であるので、第1主成分得点は麺類の消費の多さを反映すると考えられる。第2主成分ベクトルの係数は、「うどん・そば」では正、「中華麺」や「カップ麺」では負であるので、うどん・そばの消費の多さを反映すると考えられる。

pd.DataFrame(pca.components_, columns=F, index=[f'PC{i+1}' for i in range(pca.n_components_)])
生うどん・そば 乾うどん・そば パスタ 中華麺 カップ麺 即席麺 他の麺類
PC1 0.5050 0.4089 0.0169 0.5472 0.5031 -0.1329 0.0858
PC2 0.7567 0.1505 -0.0064 -0.2922 -0.5624 -0.0256 -0.0487
PC3 -0.2577 0.1220 0.0783 0.6099 -0.6228 -0.3794 0.0938
PC4 -0.3246 0.8904 -0.0693 -0.2862 -0.0487 0.0988 -0.0550
PC5 0.0214 0.0239 0.2031 0.3552 -0.1849 0.8891 -0.0829
PC6 0.0054 -0.0359 -0.6138 0.1797 -0.0219 -0.0069 -0.7675
PC7 -0.0093 0.0255 0.7555 -0.0536 0.0759 -0.1934 -0.6184

主成分分析により、各データ点を第1から第7主成分に射影したときの主成分得点を表形式で表示する。

pd.options.display.float_format = '{:.4f}'.format
pd.DataFrame(P, columns=[f'PC{i+1}' for i in range(pca.n_components_)], index=C)
PC1 PC2 PC3 PC4 PC5 PC6 PC7
北海道 153.4222 -534.3034 -152.1148 -141.8253 -420.6236 -98.0002 65.1171
青森県 1101.6032 -1496.8960 -152.2914 -252.7109 175.4282 299.5492 -12.7754
岩手県 1793.4762 -1334.3605 386.5550 -402.2430 341.3197 43.9761 -95.7998
宮城県 763.3783 -877.5282 352.0041 -76.2460 -33.1394 -91.6132 0.0315
秋田県 1172.0693 -473.2078 661.3547 767.5094 34.5863 324.9647 -63.6768
山形県 2150.6186 -95.0627 -182.2755 -10.1633 12.6520 16.8480 45.4275
福島県 783.0933 -1051.2636 -321.7877 385.6438 -435.1023 -116.8548 -166.2046
茨城県 47.2661 64.7670 221.2223 200.7490 -480.9328 -50.3860 96.3663
栃木県 642.1218 68.3046 -26.7517 -345.0909 26.5946 -115.6424 58.6243
群馬県 728.3322 594.9992 -438.8713 -701.5308 -423.2849 94.8302 114.2694
埼玉県 576.7350 361.9543 299.1710 -369.3226 -168.5329 -328.2688 33.3482
千葉県 224.5805 -90.9914 432.7337 -127.3440 98.5016 -204.3793 13.7063
東京都 -130.8999 4.4330 943.6409 62.1463 127.2547 -396.2840 27.0850
神奈川県 139.5407 275.5104 774.9649 255.8418 62.7917 -408.6177 -96.9218
新潟県 851.1234 -1217.4607 -642.1768 1196.9064 -7.8734 -103.1752 157.6930
富山県 736.6941 -612.2866 -365.4479 369.1469 -110.8816 137.8001 120.1706
石川県 547.9223 -39.0012 746.5229 -370.5084 359.9966 141.6180 27.9426
福井県 136.8892 1006.1167 418.3710 49.8097 145.8780 235.9303 68.7181
山梨県 455.0597 102.5065 123.4897 -569.4544 -249.6373 125.8385 -57.3433
長野県 549.7100 667.2875 327.6183 102.2576 -162.8597 -128.1230 -53.2867
岐阜県 -384.6726 -91.9716 -86.3750 44.1388 -158.5730 -72.7617 5.4615
静岡県 455.8524 -522.6362 874.5050 -59.6981 94.2201 -119.2878 -125.8976
愛知県 281.7116 470.1198 43.7052 -280.6772 28.3411 -151.0600 43.4152
三重県 -47.1050 439.4803 201.9089 -147.4071 67.3410 34.2168 -2.9096
滋賀県 -8.0834 411.3792 -3.9890 -344.3619 65.1042 -48.7050 -9.2709
京都府 -542.1400 1083.4309 769.0285 -147.7499 292.4058 -8.0488 -15.5865
大阪府 142.2460 -19.1879 -1296.0106 -603.7145 -160.9610 104.1794 -94.4841
兵庫県 -203.8417 31.4682 -178.0018 -176.7118 -98.2611 -95.1057 -188.2460
奈良県 -115.2237 529.8485 255.3790 462.1403 98.6741 -23.0343 -34.3392
和歌山県 -780.2281 253.3868 175.8894 435.2915 -157.0167 189.0887 -157.5332
鳥取県 -69.1043 -542.0104 -1420.7600 -331.0915 587.4943 47.9239 11.8743
島根県 -165.2852 258.9724 198.6653 72.2967 232.5081 118.8384 129.0717
岡山県 -124.8443 9.2985 -147.0644 32.0185 46.0927 -68.8385 200.0767
広島県 -588.9009 194.1220 248.4702 -254.1786 202.2912 -64.5546 78.4911
山口県 -523.1440 -6.2518 -963.2277 -90.0340 154.8302 -25.0465 -56.9032
徳島県 -623.5714 697.9746 -338.2960 61.3932 -129.3636 118.8313 18.2818
香川県 1861.8465 2898.5926 -881.5167 681.5892 130.8813 200.6197 -87.5708
愛媛県 -684.5301 866.4939 159.5305 -197.1324 -132.0062 63.2729 20.4916
高知県 -594.0772 -211.0039 -1013.0915 -374.2155 -169.9183 39.8092 22.3229
福岡県 -1183.6986 37.7734 -103.0164 -23.5416 -19.5804 -194.1989 36.1238
佐賀県 -1279.3688 -357.8023 -286.6590 443.6874 32.6912 -5.4547 -57.7901
長崎県 -1105.4465 -251.6660 274.9161 112.1158 -148.1998 175.8729 -40.7751
熊本県 -1640.3146 -571.0515 -931.9119 284.8872 329.9288 -169.0873 14.2677
大分県 -1563.1952 -255.2315 -324.8045 75.9314 57.7498 7.0017 -153.7285
宮崎県 -1431.4404 -121.9727 -200.6135 -53.1888 11.5342 33.9756 1.0191
鹿児島県 -894.9367 146.5679 130.7319 523.5690 0.6234 -11.8609 113.2129
沖縄県 -1611.2403 -701.6404 1436.6764 -168.9272 -150.9670 543.4036 48.4329

用いる主成分の数\(k\)を横軸、累積寄与率を縦軸としてグラフを描画する。第2主成分までで元データの分散の約70%、第3主成分までで元データの分散の約90%を説明できることが分かる。

fig = plt.figure(dpi=100)
ax = fig.add_subplot(1,1,1)

ax.plot(np.arange(1, pca.n_components_+1), np.cumsum(pca.explained_variance_ratio_))
ax.set_xlabel('Number of principal components (PCs)')
ax.set_ylabel('Cumulative explained variance')
fig.show()
../_images/04pca2_29_0.png

第1主成分得点を横軸、第2主成分得点を縦軸として、各都道府県庁所在地の麺類の支出金額を2次元平面上にプロットしてみる。これだけでも、地理的に近い都道府県庁所在地が近くにプロットされるのが興味深い。第1主成分得点が最も高いのは山形市(山形県)であり、周辺に他の都道府県庁所在地がプロットされていないことから、麺類の支出金額が特徴的であることが分かる。ただ、他の都道府県を全く寄せ付けない圧倒的な存在感を示すのは、第1主成分得点と第2主成分得点のいずれも高い高松市(香川県)である。

fig = plt.figure(dpi=100, figsize=(10, 10))
ax = fig.add_subplot(1,1,1)
ax.set_xlabel('PC 1')
ax.set_ylabel('PC 2')
m = ax.scatter(P[:, 0], P[:, 1], marker='.')
for i, label in enumerate(C):
    ax.text(P[i][0], P[i][1], label, ha='center', va='bottom')
ax.set_aspect('equal')
ax.grid()
fig.show()
../_images/04pca2_31_0.png

各都道府県庁所在地がどの品目に支出するのかを表すため、各品目を2,000円消費したときの点を原点からの矢印で表示した。高松市(香川県)は「生うどん・そば」に対する支出が特徴的であること、青森市(青森県)はカップ麺に対する支出が特徴的であること、山形市(山形県)は麺類全般に対して支出していることが示唆される。

fig = plt.figure(dpi=100, figsize=(10, 10))
ax = fig.add_subplot(1,1,1)
ax.set_xlabel('PC 1')
ax.set_ylabel('PC 2')
ax.set_aspect('equal')

m = ax.scatter(P[:, 0], P[:, 1], marker='.')
for i, label in enumerate(C):
    ax.text(P[i][0], P[i][1], label, ha='center', va='bottom')

for i, label in enumerate(F):
    ax.annotate(
        label, xy=(0, 0), xytext=pca.components_[:2,i] * 2000,
        color="tab:red", arrowprops=dict(arrowstyle='<-', color="tab:red")
    )
ax.grid()
fig.show()
../_images/04pca2_33_0.png

13.6. 応用例 (2): 手書き文字#

scikit-learnのsklearn.datasets.load_digitsを使い、手書き文字認識のデータセットを読み込む。

from sklearn.datasets import load_digits

X, y = load_digits(return_X_y=True)

読み込まれたデータは\(1797\)事例からなり、各事例は\(64\)\(8 \times 8\)ピクセル)次元。

X.shape
(1797, 64)
y.shape
(1797,)

データ中の先頭の事例の画像をプロットする。

fig = plt.figure(dpi=100)
ax = fig.add_subplot(1,1,1)
ax.matshow(X[0].reshape(8, -1)) 
fig.show() 
../_images/04pca2_40_0.png

入力データは\(64\)次元の特徴空間で表現されていたが、主成分分析を用いて特徴空間を\(2\)次元に圧縮する。

from sklearn.decomposition import PCA

pca = PCA(2)
P = pca.fit_transform(X)

各事例の特徴空間の次元が\(2\)次元になった。

P.shape
(1797, 2)

先ほどの画像の特徴ベクトルを\(2\)次元に圧縮したものを表示する。

P[0]
array([-1.25946627, 21.2748836 ])

各事例の特徴ベクトルが\(2\)次元に圧縮されたので、それらを\(2\)次元平面上にプロットする。そのとき、画像の数字の正解で色分けを行うと、同じ数字の画像が近くにまとまっていることが確認できる。

fig = plt.figure(dpi=100)
ax = fig.add_subplot(1,1,1)
ax.set_xlabel('PC 1')
ax.set_ylabel('PC 2')
ax.set_aspect('equal')
m = ax.scatter(P[:, 0], P[:, 1], c=y, marker='.', cmap='Set1')
fig.colorbar(m)
fig.show()
../_images/04pca2_48_0.png

第1主成分ベクトルと第2主成分ベクトルを2次元画像として可視化する。入力画像とこれらの主成分ベクトルとの内積で主成分得点(射影後の値)が計算される。

fig, axes = plt.subplots(1, 2, dpi=100)
for i, ax in enumerate(axes.flat):
    ax.matshow(pca.components_[i].reshape(8, -8))
    ax.set_title(f'PC {i+1}', y=-0.2)
../_images/04pca2_50_0.png

横軸を主成分の数、縦軸を累積寄与率としてグラフを書いてみると、おおよそ10数個の主成分ベクトルでデータ全体の分散の80%以上をカバーできることが分かる。

pca = PCA()
pca.fit(X)

fig = plt.figure(dpi=100)
ax = fig.add_subplot(1,1,1)

ax.plot(np.cumsum(pca.explained_variance_ratio_))
ax.set_xlabel('Number of principal components (PCs)')
ax.set_ylabel('Cumulative explained variance')
fig.show()
../_images/04pca2_52_0.png

13.7. 応用例 (3): 固有顔#

固有顔(Eigenface)と呼ばれる手法を体験してみよう。固有顔は顔画像のデータに対して主成分分析を行い、主成分ベクトルを求めたもので、顔認識の特徴量抽出(の次元圧縮)に用いられる。今回は、sklearn.datasets.fetch_lfw_peopleを用いて、顔画像のデータを読み込む。

from sklearn.datasets import fetch_lfw_people
faces = fetch_lfw_people(min_faces_per_person=60)

読み込まれたデータは\(1348\)件で、それぞれ、縦が\(62\)ピクセル、横が\(47\)ピクセルの画像になっている。

faces.images.shape
(1348, 62, 47)

読み込まれた画像を、その人物名と一緒に表示する。

fig, axes = plt.subplots(4, 4, dpi=100, figsize=(8, 8))
for i, ax in enumerate(axes.flat):
    ax.imshow(faces.images[i], cmap='binary_r')
    ax.set_title(faces.target_names[faces.target[i]], y=-0.2)
    ax.set_axis_off()
../_images/04pca2_58_0.png

顔画像データに主成分分析をかけて、150個の主成分ベクトルを求める。

from sklearn.decomposition import PCA

pca = PCA(150)
pca.fit(faces.data)
PCA(n_components=150)

sklearn.decomposition.PCAは、主成分分析を行う前に、データの事例ベクトルの平均を自動的に取り除く。その平均を2次元で可視化することで、先ほどのデータセットの顔画像の「平均」を見ることができる(彫りが深い…)。

fig = plt.figure(dpi=100)
ax = fig.add_subplot(1,1,1)
ax.imshow(pca.mean_.reshape(faces.images[0].shape), cmap='bone')
fig.show()
../_images/04pca2_62_0.png

第1主成分から順に、第16主成分まで可視化してみる。

fig, axes = plt.subplots(4, 4, dpi=100, figsize=(8, 8))
for i, ax in enumerate(axes.flat):
    ax.imshow(pca.components_[i].reshape(62, 47), cmap='bone')
    ax.set_title(f'PC {i+1}', y=-0.2)
    ax.set_axis_off()
../_images/04pca2_64_0.png

横軸を主成分の数、縦軸を累積寄与率としてグラフを書いてみると、だいたい30個弱の主成分ベクトルでデータ全体の分散の80%以上をカバーできることが分かる。

fig = plt.figure(dpi=100)
ax = fig.add_subplot(1,1,1)

ax.plot(np.cumsum(pca.explained_variance_ratio_))
ax.set_xlabel('Number of principal components (PCs)')
ax.set_ylabel('Cumulative explained variance')
fig.show()
../_images/04pca2_66_0.png

画像データを\(150\)次元の空間に射影してから、元の画像の空間に戻してみると、元画像を少ない情報量に圧縮した画像を得ることができる。以下の可視化では、上段が元画像、下段が圧縮後の画像である。

P = pca.transform(faces.data)
PB = pca.inverse_transform(P)

fig, axes = plt.subplots(2, 8, dpi=100, figsize=(8, 2))
for i in range(8):
    axes[0, i].imshow(faces.data[i].reshape(62, 47), cmap='binary_r')
    axes[0, i].set_axis_off()
    axes[1, i].imshow(PB[i].reshape(62, 47), cmap='binary_r')
    axes[1, i].set_axis_off()
../_images/04pca2_68_0.png