主成分分析の使いみち　次元の縮約　クラスタリング　モード分解

# はじめに
統計の教科書を読んでいるとよく出てくる主成分分析ですが、どんなところで使えばいいのかわからないのでいくつかのデータで試してみました。

私の体感ですが、統計と工学の分野で使われ方に差があると思います。

統計の分野では、3次元以上の多次元データを次元縮約して可視化する・新しい軸の解釈をする、多次元データを次元縮約してクラスタリングや機械学習の前処理に使う、重回帰分析のときに独立変数間の相関をなくす(多重共線性の解消)などです。

工学の分野(流体力学)では、例えばカルマン渦をハイスピードカメラで撮影して、その画像たちに対して主成分分析を行い、寄与率の高いモードを抜き出すことで物理現象の理解に役立てることなどに使われています。低ランク近似とも言います。
これは固有直交分解(Proper Orthogonal Decomposition, POD)と呼ばれていますが、主成分分析と同じです。

# 参考

主成分分析がわかりやすく解析されています。
[Qiita 意味がわかる主成分分析](https://qiita.com/NoriakiOshita/items/460247bb57c22973a5f0)

導出がまとめられています。この投稿では使い方のみに絞ります。
[Qiita 主成分分析について](https://qiita.com/jamojisan/items/8940cf46d4eba71b3ea3)

7科目のテストの点数を次元縮約して、どんな生徒なのか考えます。
[主成分分析 - 統計科学研究所](file:///K:/OneDrive/data/%E5%8B%89%E5%BC%B7/%E4%B8%BB%E6%88%90%E5%88%86%E5%88%86%E6%9E%90/statR_9_principal.pdf)

アヤメのデータで次元縮約してクラスタリングします。
[Qiita scikit-learnを使用した主成分分析(PCA)](https://qiita.com/kannoma/items/edff3b96b334fa19f67e)

カルマン渦の可視化画像に固有直交分解している例です。
[平邦彦, "固有直交分解による流体解析:1.基礎", ながれ, Vol.30, 2011, pp.115-123](http://www.nagare.or.jp/download/noauth.html?d=30-2rensai2.pdf&dir=48)

# 主成分分析
主成分分析は、データを軸(基底)に射影したときに分散が最も大きくなるような軸を探してくる操作です。
その目的はデータの分散共分散行列を固有値分解して、データに固有ベクトルをかけることで達成できます。
軸の直交性は保ったまま変換する直交基底変換です。
フーリエ変換なんかと同じですね。

## 使用するモジュール

```python
import numpy as np
import matplotlib.pyplot as plt
import pandas as pd
from sklearn.datasets import load_iris #あやめデータセット用
from tensorflow import keras #MINIST用
```

## 解析の流れ

1. 平均を0にする
2. 標準偏差を1にする(しないときもあります)
3. 分散共分散行列を計算する
4. 分散共分散行列を固有値分解する
5. 固有値の大きさ順に並び替える
6. データに固有ベクトルをかけて基底変換する

ここまでがどれも共通してやることです。目的に合わせて追加の操作もあります。

軸の解釈を行う場合は、元のデータと基底変換後のデータの相関を表す因子負荷量を計算します。新たしい軸が元の軸とどんな対応関係にあるかがわかります。

次元の縮約や低ランク近似で、いくつの次元を選べばいいか判断する定量的な指標として、固有値の合計に占めるn番目までの固有値の累積和を計算する累積寄与率があります。
これをもとに何次元落とすか決めます。

データを$X$とします。上の流れを式でも書きます。

```math
cov = {\rm COV}(X) \\
UWU^T = {\rm SVD}(cov) \\
{\bar X} = XU
```

$\bar X$は直交基底変換後のデータです。
$\bar X$と元のデータ$X$の相関を示す指標に因子負荷量がります。
固有値の平方根と固有ベクトルの積で求まります。
$${\sqrt W} U$$

## 実例

### 1. 多変量正規分布

手順を学ぶために多変量正規分布に主成分分析を行います。

```python
np.random.seed(0)
mean = np.array([3.0, 5.0])
cov = np.array([[1.5, 2.0], [2.0, 4.0]])
data = np.random.multivariate_normal(mean, cov, size=100)

data = data - np.mean(data,axis=0) #各次元の平均を0にする
data /= np.std(data,axis=0)
```

先に平均と標準偏差を操作してしまいましたが、これが元のデータです。
![gauss.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/455608/b4282a5e-a326-9f51-23bb-0163b3b965a0.png)

```python
cov = np.cov(data, rowvar=False) #分散共分散行列
u,w,_ = np.linalg.svd(cov) #固有値分解
pca_data = data.dot(u)

pca_ax = u.dot(u) #赤矢印

fig,axes = plt.subplots(figsize=(6,8))
axes.scatter(pca_data.T[0],pca_data.T[1]) #元データ回転
axes.scatter(pca_data[:,0].T,-2.0*np.ones(pca_data.shape[0]),alpha=0.3) #PC1射影
axes.scatter(3.8*np.ones(pca_all.shape[0]),pca_data[:,1].T,alpha=0.3)
axes.set_aspect("equal")
axes.set_ylim(-4.0,4.0)
axes.set_xlim(-4.0,4.0)
axes.quiver(0,0,pca_ax[0,0],pca_ax[0,1],color = "red",angles = "xy", scale_units = "xy",scale = 0.5)
axes.quiver(0,0,pca_ax[1,0],pca_ax[1,1],color = "red",angles = "xy", scale_units = "xy",scale = 0.5)
```
上の図の赤矢印は$U$をプロットしたものです。

直交基底変換後
![gauss_kaiten_touei.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/455608/cbb40600-4a6e-7442-a180-acc8bf4b5919.png)

モード分解はこのように次元を落として逆変換します。
フーリエ変換して、周波数空間でフィルタをかけて、逆フーリエ変換してもとに戻すのと同じ感じです。

```python
pca_approx = np.copy(pca_data)
pca_approx[:,1:] = 0.0 #低ランク近似
recon = pca_approx.dot(u.T) #逆変換

fig,axes = plt.subplots(figsize=(6,8))
axes.scatter(data.T[0],data.T[1])
axes.scatter(recon.T[0],recon.T[1])
axes.quiver(0,0,u[0,0],u[0,1],color = "red",angles = "xy", scale_units = "xy",scale = 0.5)
axes.quiver(0,0,u[1,0],u[1,1],color = "red",angles = "xy", scale_units = "xy",scale = 0.5)
axes.set_aspect("equal")
axes.set_ylim(-4.0,4.0)
axes.set_xlim(-4.0,4.0)
```

![gauss_kaiten_touei_saikousei.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/455608/9ec2b032-660d-ba01-a126-51ace24ab347.png)

累積寄与率も見ます。
次元の縮約をするときに、いくつの次元を残せばだいたいデータを説明できるかを定量的に知るための指標です。

```python
np.sum(w[0])/np.sum(w)

>>> 0.9080218762749276
```

### 2. 次元の縮約　テストの点数

統計科学研究所のテストデータを使用します。9科目166人分のデータです。
9次元あるので、生徒がどんな傾向にあるのか判断するのが難しいです。
そこで2次元まで落として、軸を解釈して生徒の傾向を掴みます。
[統計科学研究所 統計データ](https://statistics.co.jp/reference/statistical_data/statistical_data.htm)

#### 2-1. データの読み込み

```python
csv_data = pd.read_csv("seiseki.csv")

csv_data
>>>
kokugo	shakai	sugaku	rika	ongaku	bijutu	taiiku	gika	eigo
0	30	43	51	63	60	66	37	44	20
1	39	21	49	56	70	72	56	63	16
2	29	30	23	57	69	76	33	54	6
3	95	87	77	100	77	82	78	96	87
4	70	71	78	67	72	82	46	63	44
...	...	...	...	...	...	...	...	...	...
161	82	78	80	88	80	69	83	78	90
162	0	8	2	9	5	18	42	2	1
163	45	26	29	24	31	57	68	40	27
164	73	31	43	32	59	64	82	48	56
165	60	85	89	80	56	79	81	45	85

166 rows × 9 columns
```

#### 2-2. 散布図行列・相関行列

データをいきなり主成分分析にかけてしまう前に、散布図行列を見てどんなデータなのか確認します。

```python
pd.plotting.scatter_matrix(csv_data, figsize=(15, 15))
```
![散布図行列.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/455608/9b58b301-5e03-1097-36dd-b6c88e0c799c.png)

正直これを見ただけでも主成分分析の結果がだいたいどうなるかわかってしまいますね。
国語、社会、理科、音楽、美術、技術家庭科、英語はそこそこ相関があります。
体育だけ他の教科と相関が低いです。

相関行列も見ます。

```python
csv_data.corr().style.background_gradient(axis=None, cmap="jet")
```
![相関行列.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/455608/50fb8f57-8695-8143-88fc-03a666d7aabf.png)

#### 2-3. 主成分分析

主成分分析を行います。

```python
col = csv_data.columns #列名
data = csv_data.to_numpy().astype("float64")

data -= np.mean(data,axis=0) #各列の平均を0にする
data /= np.std(data,axis=0) #各列を標準偏差で正規化

cov = np.cov(data, rowvar=False) #分散共分散行列
u,w,_ = np.linalg.svd(cov) #固有値分解 u:固有値ベクトル w:固有値

sum_w = np.sum(w)
cr = []
for w_ in w:
    cr.append(w_/sum_w)
ccr = np.cumsum(cr) #累積寄与率
```

まずは累積寄与率を見ます。

```python
fig,axes = plt.subplots()
axes.plot(list(range(1,len(ccr)+1)),ccr,marker="o")
axes.set_xlabel("Components number")
axes.set_ylabel("Cumulative contribution ratio")
axes.set_xlim(0,10)
```

![seiseki_ruisekikiyo.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/455608/802e8c1b-cad3-cb3c-2ee5-09e810fe907d.png)

2次元まで落としても、80%くらいの情報は残っているようです。

主成分得点を見ます。

```python
pca = data.dot(u)

fig,axes = plt.subplots(figsize=(6,8))
axes.scatter(pca.T[0],pca.T[1])
axes.set_aspect("equal")
axes.set_xlim(-6.0,6.0)
axes.set_ylim(-6.0,6.0)
axes.set_xlabel("PC1")
axes.set_ylabel("PC2")
```

![seiseki_biplot.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/455608/74d0759f-3e9c-9e36-0dce-04aec200f550.png)

9次元のデータが2次元になりました。しかし、PC1とPC2の意味がわかりません。
そこで元のデータと変換後のデータの相関を表す因子負荷量を計算します。

```python
fl = np.sqrt(w) * u
```

```python
#FC1
fig,axes = plt.subplots()
axes.scatter(list(range(fl.shape[0])),fl[:,0])
plt.xticks(list(range(fl.shape[0])),list(col))
axes.set_ylabel("FC1")
plt.xticks(rotation = 20)
```

第一主成分軸と関係が深い教科の値が大きくでます。
体育以外の教科と関連が深いようです。
![seiseki_fc1_1d.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/455608/51fde881-1073-5ea4-18e6-902f3ea7dd67.png)

```python
#FC2
fig,axes = plt.subplots()
axes.scatter(list(range(fl.shape[0])),fl[:,1])
plt.xticks(list(range(fl.shape[0])),list(col))
plt.xticks(rotation = 20)
axes.set_ylabel("FC2")
```
![seiseki_fc2_1d.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/455608/cd4821ac-6b05-ae5c-aa57-3a80910a2531.png)

符号は反転ですが、値の絶対値の大きさに意味があります。
第二主成分軸は体育が強く関係します。

同時にプロットしてもわかりやすいです。

```python
def circle(n,r=1):
    """
    円を描く nは分割数
    """
    rad = np.linspace(0.0,2.0*np.pi,n)[1:] #0と2piが重ならないように
    return np.cos(rad),np.sin(rad)

c = circle(100)

fig,axes = plt.subplots(figsize=(6,8))
axes.scatter(c[0],c[1],marker=".",color="k",alpha=0.2) #円の描画
for i in range(u.shape[0]):
    axes.scatter(fl[i,0],fl[i,1],label=col[i])

axes.set_aspect("equal")
axes.legend(loc="best")
axes.set_xlabel("FC1")
axes.set_ylabel("FC2")
```
![seiseki_fc_2d.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/455608/bac958a3-bb79-99f7-b532-f670fc289801.png)

体育以外の科目は第一主成分軸と関係が深い一方、第二主成分軸にはあまり寄与しないようです。
体育はその反対です。

ここでもう一度主成分得点のプロットを見ます。
右に行くほど筆記試験(体育以外)の得点が高く、下に行くほど体育が得意な生徒ということになります。
このように次元を縮約して、新しい軸に意味付けすることでデータを解釈しやすくなります。
![seiseki_biplot.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/455608/040ed9b8-17e2-ed81-ec01-fdee297ede6b.png)

### 3. クラスタリング　あやめのデータ
次は3種類(セトサ種, バージカラー種, バージニカ種)のあやめの4つの特徴(ガクの長さ, ガクの幅, 花弁の長さ, 花弁の幅)を記録したデータでクラスタリングを行います。

#### 3-1. データの読み込み
データはsklearnから持ってきます。

```python
iris = load_iris()

iris.keys()
>>> dict_keys(['data', 'target', 'target_names', 'DESCR', 'feature_names', 'filename'])

iris["target_names"]
>>> array(['setosa', 'versicolor', 'virginica'], dtype='<U10')

iris["feature_names"]
>>> ['sepal length (cm)',
     'sepal width (cm)',
     'petal length (cm)',
     'petal width (cm)']

iris_df = pd.DataFrame(iris["data"],columns=iris["feature_names"])

iris_df
>>>
sepal length (cm)	sepal width (cm)	petal length (cm)	petal width (cm)
0	5.1	3.5	1.4	0.2
1	4.9	3.0	1.4	0.2
2	4.7	3.2	1.3	0.2
3	4.6	3.1	1.5	0.2
4	5.0	3.6	1.4	0.2
...	...	...	...	...
145	6.7	3.0	5.2	2.3
146	6.3	2.5	5.0	1.9
147	6.5	3.0	5.2	2.0
148	6.2	3.4	5.4	2.3
149	5.9	3.0	5.1	1.8
150 rows × 4 columns
```

#### 3-2. 散布図行列・相関行列
ここでも散布図行列と相関行列を見ます。

```python
pd.plotting.scatter_matrix(iris_df, figsize=(15, 15),c=iris["target"],marker="o")
```
![Iris_分散共分散行列.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/455608/6530d69d-22be-cc17-606d-b2bab5371209.png)

なんとなく種ごとにプロットがかたまっていますね。

```python
iris_df.corr().style.background_gradient(axis=None, cmap="jet")
```
![相関行列.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/455608/963b3652-fce7-57d5-22fa-5042080745c3.png)

#### 3-3. 主成分分析

データを正規化して主成分分析をやっていきます。

```python
data = iris["data"]
data -= np.mean(data,axis=0) #各列の平均を0にする
data /= np.std(data,axis=0) #各列を標準偏差で正規化

cov = np.cov(data, rowvar=False) #分散共分散行列
u,w,_ = np.linalg.svd(cov) #固有値分解 u:固有値ベクトル w:固有値

sum_w = np.sum(w)
cr = []
for w_ in w:
    cr.append(w_/sum_w)
ccr = np.cumsum(cr) #累積寄与率

fl = np.sqrt(w) * u #因子負荷量
```

グラフの描画は先程と同様なので省略します。

累積寄与率を見ると、2次元だけで95%以上の情報が表現できます。

![iris_kiyoritu.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/455608/b97a0173-c2a9-7491-c5e1-fc7a50b4ea1a.png)

因子負荷量を見ます。
sepal length(ガクの長さ)とpetal length(花弁の長さ), petal width(ガクの幅)は第一主成分軸に強い相関があります。
sepal width(ガクの幅)は第二主成分軸と強い相関があります。

![fc_2d.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/455608/151b0a96-325f-1406-5841-dbf407720ea7.png)

4次元を2次元に落としました。種別にまとまっています。
今回は答えがわかっている状態ですが、種のデータがなく特徴だけのデータでも、k平均法なり線形判別分析なんなりでクラスタリングができます。

![pca.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/455608/be074ae2-15ff-6aff-5464-7dba6ba34bc3.png)

### 4. モード抽出・低ランク近似　MINIST

今度は手書き文字データ(MINIST)の数字の「3」についてモード抽出を行います。
MINISTは28x28のデータなので784次元のデータと見なせます。
これを例えば2次元などに落として逆変換することで、特徴的なモードを見ることができます。

#### 4-1. データの読み込み
kerasからデータを持ってきます。

```python
mnist = keras.datasets.mnist
(train_images, train_labels), (test_images, test_labels) = mnist.load_data()

num_3 = train_images[train_labels==3][0:1000] #3だけ

flatten_imgs = num_3.reshape(1000,28*28) #28x28の画像を1次元配列にする
````

#### 4-2. 主成分分析　モード抽出・低ランク近似

```python
flatten_imgs = flatten_imgs - np.mean(flatten_imgs,axis=0) #各列の平均を0にする

cov = np.cov(flatten_imgs, rowvar=False) #分散共分散行列
u,w,_ = np.linalg.svd(cov) #固有値分解 u:固有値ベクトル w:固有値

sum_w = np.sum(w)
cr = []
for w_ in w:
    cr.append(w_/sum_w)
ccr = np.cumsum(cr) #累積寄与率
```

784次元までありますが、100次元までを表示しています。
100次元もあればほぼすべての情報が表現できるようです。
今回は試しに3次元と100次元で近似してみます。

![寄与率.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/455608/d586f19c-e20d-04fc-4c1a-b56024b713ba.png)

```python
pca = flatten_imgs.dot(u)

rank = 3
pca_approx = np.copy(pca)
pca_approx[:,rank:] = 0.0

rebuilt = pca_approx.dot(u.T)
rebuilt_reshape = rebuilt.reshape(1000,28,28)
```

数枚の画像を表示します。
上列がランク3近似、下列がランク100近似です。

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/455608/b46a24ce-c92f-fa82-ee83-12d2f411e5e9.png)

ランク3だとデータ間の差が少なく、特徴的なモードのみが残っています。
ランク100だとデータ間の差がはっきりわかり、元のデータがよく再現されています。

# まとめ
主成分分析の使いみちを3つ試してみました。
次元の縮約、クラスタリング、モード分解のどれも同じような手順で行うことができます。