2. Pythonで綴る多変量解析　3-1. 主成分分析(scikit-learn)

**主成分分析**とは、多数の変量を、少数の新変量に要約する技法です。すなわち、いくつもの変数がもっている情報を圧縮して新たな変数を合成し、それでデータを再構成することによって全体の傾向や個体の特徴を顕在化させるテクニックといえます。この新たな変数のことを**主成分**と呼びます。

かりに、数学・理科・社会・英語・国語という5教科のテストの得点が1クラス20名分あったとします。そこで、どのくらい学力があるのか知りたいとき、普通は「合計点」でみるということをやると思います。5教科を変数$x_{1}, x_{2}, x_{3}, x_{4}, x_{5}$として式に表わしてみます。

```math
合計点 = 1×x_{1} + 1×x_{2} + 1×x_{3} + 1×x_{4} + 1×x_{5}
```

わざわざ $1×$ と書いていますが、5教科にそれぞれ同じ重み$1$をかけて足し上げたのが合計点です。ここで何が起こっているのかというと、もともと5教科＝5次元あった情報を「学力」という1次元に要約したことになります。この合計点で学力をみるというやり方は、その前提として学力のある人はどの教科も得点が良いだろうという思いがあります。これを言い換えると、各教科の間には相関があるだろう、となります。

そこで、具体的に1クラス20名分のダミーデータを用意して、ひとまずデータ全体の様子を見てみましょう。

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

```python
# 数値計算やデータフレーム操作のためのライブラリ
import numpy as np
import pandas as pd

# グラフを描画するためのライブラリ
import matplotlib.pyplot as plt
%matplotlib inline

# 機械学習ライブラリ
import sklearn
from sklearn.decomposition import PCA
```

####**⑵ データを用意する**

```python
# 20人分のテスト得点データ
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 = ['数学', '理科', '社会', '英語', '国語'])
```

![002_0301_002.PNG](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/502109/1cc152d4-2aa4-2685-b2ae-18764c92db3e.png)

####**⑶ 散布図を描いてデータを概観する**

```python
# PandasのPlotメソッドをインポート
from pandas import plotting

# 散布図を描画
plotting.scatter_matrix(df, alpha=0.5, figsize=(8, 8))
plt.show()
```

散布図を描くコード`pandas.plotting.scatter_matrix(frame, alpha, figsize)`のうち、引数の`alpha`は描画色の透過率（0～1）、`figsize`は描画サイズ（幅, 高さ）をインチ単位で指定しています。

![002_0301_003.PNG](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/502109/476d4305-28b0-5498-bf74-13ce77d4a30f.png)

いかがでしょう。たとえば理科と数学など、右肩上がりの直線的な分布がみてとれます。理科の得点が高いと数学の得点も高い傾向がある、つまりデータの間に相関がある。とすれば、情報をもっとコンパクトに圧縮できるのではないか、というのが**主成分分析**の考え方です。

主成分分析 (PCA)
---
主成分分析の狙いは、情報をなるべく失わずにデータを圧縮するということです。次元数を減らすことで生じる情報の損失を最小限におさえたい、言い換えれば元の情報量を最大限に残してくれる新たな変数を求めます。
では実際に主成分分析はこんなステップを踏んでやっていきますというところを確認していきます。

まず、5教科＝5変数のベクトルを次のようにおきます。

```math
  x = \left(
    \begin{array}{ccc}
      x_{1} \\
      x_{2} \\
      x_{3} \\
      x_{4} \\
      x_{5} 
    \end{array}
  \right)
```

**第1主成分の計算**
主成分を、次のように「各変数に係数をつけた各成分の値を足してやったもの」と定義します。

```math
z_{1} = w_{11}x_{1} + w_{12}x_{2} + w_{13}x_{3} + w_{14}x_{4} + w_{15}x_{5} = w_{1}･x
```

右辺の$w_{1}$をベクトルに表わすと次のようになりますが、先ほどの合計点の例ではここに全て$1$が入っていました。

```math
  w_{1} = \left(
    \begin{array}{ccc}
      w_{11} \\
      w_{12} \\
      w_{13} \\
      w_{14} \\
      w_{15} 
    \end{array}
  \right)
```

では、どのような$w_{1}$を探せばいいかといえば、くり返しになりますが「情報量が最大になる$w_{1}$」です。
この情報量とは何なのか、主成分分析では「情報＝分散」ということになっています。分散とはデータの散らばり具合いですが、ではなぜ散らばりが大きいと情報量が大きいといえるのか。あくまでたとえ話ですが、授業の終りにやるごく簡単な小テストで全員が10点満点だったとします。つまり分散は$0$という状態ですが、こうなると個体を特徴づける情報が全くないことになります。対して、例題の5教科のように分散がそれなりにあれば、たとえば90点以上なら良くできたとか、20点台となると良くないとか評価できるわけです。そういうわけで「分散が最大＝情報が最大」ということができます。

ただし、分散が最大といっても、実はいくらでも大きくできるのです。
$w_{1}$から$w_{5}$まで全部の値を100倍すると、$z_{1}$も100倍になって、分散は10000倍になります。際限なく大きくできてしまう、そういうことではなくて、知りたいのは教科ごとの重みの割合です。どういうふうな比率で割りふると、最も情報量が大きくなるのかを知りたいわけです。だから大きさを一定にしておきます。

```math
\|w_{1}\| = 1
```
重みを全部足すと$1$になる、という一定のルールのもとに情報量$V[z_{1}]$が最大になる$w_{1}$を探します。つまり、どちらの方向にどういう割合で足すと情報量が最大になるのかが見つかると、これを**第1主成分**といいます。
これをくり返していくのです。

**第2主成分の計算**
ただし、もう情報量が最大なのは$w_{1}$ということはわかっているので、今度は$w_{1}$とは違うタイプの重みのつけ方で情報量が最大になるものを探したい。そこで、$w_{2}$は$w_{1}$と違う方向を向いていて欲しいので、条件を1つ加えます。

```math
\|w_{2}\| = 1, w_{2}\perp{w_{1}}
```
$w_{2}$は$w_{1}$と直交する、という条件です。これで$w_{2}$は$w_{1}$とは別の種類の情報をもつことになります。
次いで$w_{3}$のときは、$w_{1}$とも$w_{2}$とも直交するという条件をつけ、さらに$w_{4}$のときは･･･というふうに条件を加えながらくり返していくことによって、次のようなものが得られます。

```math
  \left(
    \begin{array}{ccc}
      z_{1} \\
      \vdots \\
      z_{5} 
    \end{array}
  \right)
  = \left(
    \begin{array}{cccc}
      w_{11} & \ldots & w_{15} \\
      \vdots & \ddots & \vdots \\
      w_{51} & \ldots & w_{55}
    \end{array}
  \right)x
```

```math
z=Wx
```
こういう数式をつかって、数万とか数十万次元のデータから、もとの情報がなるべく残るようにして数百次元くらいまでに圧縮するというのが**主成分分析**です。つまり、第k主成分は、第k番目にデータのばらつきが大きい方向というわけです。
ちなみに、各変数の重み$w_{1}, w_{2}, w_{3}, w_{4}, w_{5}$のことを**主成分負荷量**とも呼びます。
ということで、いよいよ機械学習ライブラリのscikit-learnを利用して、この**主成分$z_{1}$**と**主成分負荷量$w_{1}, w_{2}, w_{3}, w_{4}, w_{5}$**を求めます。


####**⑷ 主成分分析を実行する**

```python
# モデルのインスタンスを作成
pca = PCA()

# データをもとにモデルを作成
pca.fit(df)

# モデルにデータを当てはめる
values = pca.transform(df)
```

主成分分析は英語でPrincipal Component Analysisというので、PCAと略して呼ばれます。
まず、➀モデルのひな型となるインスタンスを作成して、➁そのインスタンスに`fit`関数でデータを渡してやるとモデルが生成されます。➂このモデルに改めてデータを当てはめてやると主成分ごとの得点が算出される、という3ステップです。

![002_0301_004.PNG](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/502109/69b9d8eb-cfbd-0a45-75f8-3a432e5e6d87.png)

なんとも見づらいので、データフレームに変換します。

```python
df_pca = pd.DataFrame(data = values,
                      columns = ["主成分{}".format(x+1) for x in range(len(df.columns))])
```

![002_0301_005.PNG](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/502109/4346dc1a-ca5a-2b86-755e-1128e9377d17.png)

第1主成分から第5主成分まであって、各生徒が5つずつ得点をもっています。この得点を**主成分得点**といいます。もとが5教科で5次元ですから、主成分もまた最も多くて5次元になります。

####**⑸ 寄与率を計算する**

```python
# モデルpcaの寄与率を算出
ev_ratio = pca.explained_variance_ratio_

# 寄与率をデータフレームに変換
df_evr = pd.DataFrame(data = ev_ratio, 
                      columns = ['寄与率'],
                      index = ["主成分{}".format(x+1) for x in range(len(df.columns))])
```

![002_0301_006.PNG](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/502109/561c9b1e-022b-c83e-98c9-c1f62732d46d.png)

寄与率とは、主成分それぞれの説明力を表わす指標です。データがもともと持っている情報量＝分散に対して、その主成分の分散＝情報量が何％にあたるかという割合なので$0≦c≦1$の値をとります。とにかく第1主成分が情報量最大で、それ以外で情報量最大なのは、そのまた次に最大なのは･･･というふうに探し出していくので、寄与率は第1主成分が最も大きく順に小さくなり、すべての寄与率の合計は$1$になります。これを**累積寄与率**というグラフにしてみるとよくわかります。

```python
# 寄与率を累積する
cc_ratio = np.cumsum(ev_ratio)

# 0を連結
cc_ratio = np.hstack([0, cc_ratio])

# グラフを描画
plt.plot(cc_ratio, "-o")
plt.xlabel("主成分")
plt.ylabel("累積寄与率")
plt.grid()
plt.show()
```

第1主成分だけで寄与率が90%を超える結果となっているので、第2主成分以降はいらない気もしますが、必ずしも常にこういうふうになるわけではありません。
![002_0301_007.PNG](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/502109/1ca8be35-f534-f3c1-5f04-bd3d2ba98c3a.png)

####**⑹ 主成分を解釈する**

```python
# 主成分の分散を算出
eigen_value = pca.explained_variance_
pd.DataFrame(eigen_value, 
             columns = ["分散"], 
             index = ["主成分{}".format(x+1) for x in range(len(df.columns))])
```

![002_0301_008.PNG](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/502109/1c42fbc0-4532-aef2-a49a-365d7a855a10.png)

分散の大きさが寄与率を反映していることがわかります。第1主成分の分散は、他の主成分に対して圧倒的に大きく、その分情報量もたくさん持っています。

```python
# 主成分負荷量を算出
eigen_vector = pca.components_
# データフレームに変換
pd.DataFrame(eigen_vector,
             columns = [df.columns],
             index = ["主成分{}".format(x+1) for x in range(len(df.columns))])
```

![002_0301_009.PNG](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/502109/3087a981-151f-e8cd-15b3-27abb0813820.png)

主成分にはどんな意味があるのでしょうか。
**第1主成分**：5教科全部にマイナスの符号がついています。これは合計点を逆転したようなもので、要は全部高いか全部低いかという方向にデータが一番ばらついているということです。マイナスがついているので、合計点が高ければ高いほど主成分得点は小さくなります。その中でも特に英語の係数が大きくなっていて、つまり英語の得点が少し違うだけで主成分得点は大きく違ってきます。
**第2主成分**：英語・国語がマイナスで大きく、数学・理科がプラスで大きくなっています。数学・理科の得点が高いと主成分得点は大きくなり、また英語・国語の得点が低いとやはり主成分得点は大きくなります。いわば理系か文系かという方向にデータがばらけており、理系科目は得意で文系科目は苦手という人ほど主成分得点は大きくなります。
**第3主成分**：国語が飛び抜けていてマイナスになっています。つまり国語ができるか否かという方向にデータがばらけており、国語の得点が高ければ高いほど主成分得点は小さくなります。

まとめると、一番データがばらついている方向（**第1主成分**）は合計点が高いか低いか。その次（**第2主成分**）は文系寄りか理系寄りか、更にその次にばらつきが大きい方向（**第3主成分**）というのが国語ができるか否か、そんな感じにデータがばらけているようです。
この3つの主成分をつかえば、3次元に圧縮する中では一番豊かに情報が手に入るということが、この分析からわかります。

####**⑺ 合計点と第1主成分の得点を比較してみる**
終りに、単純に5教科の点数を足し上げた合計点と、5教科分の情報量の9割方は保有している第1主成分の得点とで、生徒の学力ランキングがどう違ってくるのか確かめておきたいと思います。

```python
## 合計点によるランキングを作成
# 合計点を計算
sum = np.sum(np.array(df), axis=1)
# 20×1の2次元配列に変換
sum.reshape(len(sum), 1)
# データフレームに変換
df_sum = pd.DataFrame(sum,
                      columns = ['合計点'],
                      index = ["ID{}".format(x+1) for x in range(len(df.index))])
# 降順ソート
df_sum_rank = df_sum.sort_values('合計点', ascending=False)

## 第1主成分得点によるランキングを作成
# 主成分1を抽出
df_PC1 = df_pca["主成分1"]
# arrayに変換
pc1 = np.array(df_PC1)
# 20×1の2次元配列に変換
pc1 = pc1.reshape(len(pc1), 1)
# IDをつける
df_pca = pd.DataFrame(pc1,
                      columns = ['主成分1'],
                      index = ["ID{}".format(x+1) for x in range(len(df.index))])
# 昇順ソート
df_pca_rank = df_pca.sort_values('主成分1')
```
![002_0301_012.PNG](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/502109/b3ad36c1-5612-6ade-315e-3011e9f66315.png)

**第1主成分**は合計点が高いか低いかということですから、上位グループと下位グループの順位は一致しています。そのあいだの中間層は多少前後していますが、おおむね「学力のある人はどんな教科もできるんじゃないかな」という合計点の気持ちは適っていることがわかりました。

---

さて次節は、scikit-learnを使わずに、さらに主成分分析の計算のしくみをひも解いてみたいと思います。