Python 機械学習プログラミングの備忘録その５（主成分分析）

#1.はじめに
**Python機械学習プログラミング[第2版]** という本を買って勉強を始めました。効率的な学習を進める為には、勉強した内容をアウトプットするのが効果的ということで、勉強した内容を備忘録の形でここに残すことを進めています。

なお、本ではコードを一部 jupyter notebook形式で記載しているようですが、ここでは.py形式でも実行できるような形式で記載したいと思います。

今回は、第５章次元削減の中の主成分分析です。

#2.主成分分析（wineデータ）
主成分分析の目的は、分散が最大となる方向を複数見つけ出し、寄与率の高いものに絞り込み（次元削減し）、データを新たな空間に射影することです。

まず、テキストにあるコードを実行してみましょう。

```python
import pandas as pd
from sklearn.model_selection import train_test_split
from sklearn.preprocessing import StandardScaler
import numpy as np
import matplotlib.pyplot as plt

# データセットの読み込みとサマリー表示
df_wine = pd.read_csv('https://archive.ics.uci.edu/ml/'
                      'machine-learning-databases/wine/wine.data',
                      header=None)

df_wine.columns = ['Class label', 'Alcohol', 'Malic acid', 'Ash',
                   'Alcalinity of ash', 'Magnesium', 'Total phenols',
                   'Flavanoids', 'Nonflavanoid phenols', 'Proanthocyanins',
                   'Color intensity', 'Hue',
                   'OD280/OD315 of diluted wines', 'Proline']
print(df_wine.head())
print()

# データセットのトレーニングデータとテストデータへの分割(7:3)
X, y = df_wine.iloc[:, 1:].values, df_wine.iloc[:, 0].values
X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.3, 
                     stratify=y, random_state=0)  
print('X_train.shape = ', X_train.shape)
print('y_train.shape = ', y_train.shape)

# Xデータの標準化
sc = StandardScaler()
X_train_std = sc.fit_transform(X_train)
X_test_std = sc.transform(X_test)

# 固有値と固有ベクトルの計算
cov_mat = np.cov(X_train_std.T)  # 共分散行列を計算
eigen_vals, eigen_vecs = np.linalg.eig(cov_mat)  # 共分散行列から固有値と固有ベクトルを計算
print('\nEigenvalues \n%s' % eigen_vals)  # 固有値の表示
print('\nEigenvectors \n%s' % eigen_vecs)  # 固有ベクトルの表示
```
![スクリーンショット 2019-09-29 11.05.18.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/209705/46ba4b50-667b-de7d-5c5f-350f7008332d.png)
![スクリーンショット 2019-09-29 11.06.00.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/209705/d45db607-a284-5a59-72a2-13ce3797f1f0.png)
wineデータセットを読み込み、トレインとテストに7:3で分割します。ここで実際に使うのはトレインデータのみです。X_train.shape = (124, 13)なので、13次元のデータが124個あるわけです。

X_trainデータを標準化し、共分散行列(cov_mat)から固有値(Eigenvalues)と固有ベクトル(Eigenvectors)を計算します。13次元のデータなので、固有値は13個、それに対応する固有ベクトルは、13次元のものが13個得られます。

```python
# 主成分のパレート図表示
tot = sum(eigen_vals)
var_exp = [(i / tot) for i in sorted(eigen_vals, reverse=True)]
cum_var_exp = np.cumsum(var_exp)
plt.bar(range(1, 14), var_exp, alpha=0.5, align='center',
        label='individual explained variance')
plt.step(range(1, 14), cum_var_exp, where='mid',
         label='cumulative explained variance')
plt.ylabel('Explained variance ratio')
plt.xlabel('Principal component index')
plt.legend(loc='best')
plt.tight_layout()
plt.show()

# 射影行列wの作成
eigen_pairs = [(np.abs(eigen_vals[i]), eigen_vecs[:, i])  # 固有対（固有値と固有ベクトルのセット）を作成
               for i in range(len(eigen_vals))]  
eigen_pairs.sort(key=lambda k: k[0], reverse=True)  # 固有対のソート
w = np.hstack((eigen_pairs[0][1][:, np.newaxis],
               eigen_pairs[1][1][:, np.newaxis]))  # 上位2つのベクトルを取得
print('Matrix W:\n', w)  # 射影行列w(13行×2列)を表示

# データセットを射影行列wで変換
X_train_pca = X_train_std.dot(w)

# 変換したデータをプロット
colors = ['r', 'b', 'g']
markers = ['s', 'x', 'o']

for l, c, m in zip(np.unique(y_train), colors, markers):
    plt.scatter(X_train_pca[y_train == l, 0], 
                X_train_pca[y_train == l, 1], 
                c=c, label=l, marker=m)
plt.xlabel('PC 1')
plt.ylabel('PC 2')
plt.legend(loc='lower left')
plt.tight_layout()
plt.show()
```
![スクリーンショット 2019-09-29 11.22.09.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/209705/70fd83b7-a8b7-bb2a-b79f-75563891959a.png)
13個ある固有値の個々の大きさは、全体への寄与度を表すので、パレート図にしてみます。折れ線グラフは累積寄与度を表すので、上位2つの固有値で、全体の6割近くに寄与することが分かります。

固有対（固有値と固有ベクトルをセットしたもの）をソートし、上位２つの固有ベクトルを抽出し、射影行列W（13行×2列)を作成します。

13次元のトレイニングデータと2次元の射影行列の内積を取ることで、トレイニングデータが2次元に変換されます。これで、2次元グラフにプロット出来るようになるわけです。

実際に、2次元グラフにプロットすると、一目でこれなら線形分割できそうな感じが分かると思います。

#3.主成分分析（irisデータ）
テキストには載っていませんが、シンプルなirisデータセットでも、やってみます。せっかくなので、共分散行列と固有値、固有ベクトルの関係も確認してみましょう。

```python
import pandas as pd
from sklearn.model_selection import train_test_split
from sklearn.preprocessing import StandardScaler
import numpy as np
import matplotlib.pyplot as plt
from sklearn import datasets

# irisデータセットの読み込み
iris = datasets.load_iris()
X = iris.data
y = iris.target 

X_train, X_test, y_train, y_test =     train_test_split(X, y, test_size=0.3, 
                     stratify=y,
                     random_state=0)  
print('X_train.shape = ', X_train.shape)
print('y_train.shape = ', y_train.shape)

# Xデータの標準化
sc = StandardScaler()
X_train_std = sc.fit_transform(X_train)
X_test_std = sc.transform(X_test)

# 固有値と固有ベクトルの計算
cov_mat = np.cov(X_train_std.T)  # 共分散行列を計算
eigen_vals, eigen_vecs = np.linalg.eig(cov_mat)  # 共分散行列から固有値と固有ベクトルを計算
print('\ncov_mat \n%s' % cov_mat)  # 共分散行列を表示
print('\nEigenvalues \n%s' % eigen_vals)  # 固有値の表示
print('\nEigenvectors \n%s' % eigen_vecs)  # 固有ベクトルの表示

array = np.array([[eigen_vals[0], 0, 0, 0],
                  [0, eigen_vals[1], 0, 0],
                  [0, 0, eigen_vals[2], 0],
                  [0, 0, 0, eigen_vals[3]]])  # eigen_valsからarrayを作成 
print('\narray \n%s' % array) 
print('\nEigenvectors(-1) \n%s' % np.linalg.inv(eigen_vecs))  # eigen_vecsの逆行列
results = eigen_vecs.dot(array).dot(np.linalg.inv(eigen_vecs))
print('\nEigenvectors * array * Eigenvectors(-1) = \n%s' % results)
```
![スクリーンショット 2019-09-30 10.10.11.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/209705/c01ca5c2-47a4-3728-6d19-5efeac712a8a.png)
train.shape = (105, 4)なので、トレインデータは4次元です。なので、固有値（Eigenvalues)は4個、固有ベクトル(Eigenvectors)は4次元が4個です。

さて、共分散行列(cov_mat)、固有値(Eigenvalues)、固有ベクトル(Eigenvectors)の関係を見てみましょう。

固有値(Eigenvalues)を対角に入れ、その他を0とした行列を array、固有ベクトル(Eigenvectors)の逆行列を Eigenvectors(-1)とします。

そうすると、**共分散行列 = 固有ベクトル ＊ array ＊ 固有ベクトルの逆行列**　という関係にあることが分かると思います （＊は内積）。ちなみに、ある正方行列をこの様な形に分解することを**固有値分解**と言います。


さて、この後は、

```python
# 主成分のパレート図表示
tot = sum(eigen_vals)
var_exp = [(i / tot) for i in sorted(eigen_vals, reverse=True)]
cum_var_exp = np.cumsum(var_exp)

plt.bar(range(1, 5), var_exp, alpha=0.5, align='center',
        label='individual explained variance')
plt.step(range(1, 5), cum_var_exp, where='mid',
         label='cumulative explained variance')
plt.ylabel('Explained variance ratio')
plt.xlabel('Principal component index')
plt.legend(loc='best')
plt.tight_layout()
plt.show()

# 射影行列wの作成
eigen_pairs = [(np.abs(eigen_vals[i]), eigen_vecs[:, i])  # 固有対（固有値と固有ベクトルのセット）を作成
               for i in range(len(eigen_vals))]  
eigen_pairs.sort(key=lambda k: k[0], reverse=True)  # 固有対のソート
w = np.hstack((eigen_pairs[0][1][:, np.newaxis],
               eigen_pairs[1][1][:, np.newaxis]))  # 上位2つのベクトルを取得
print('Matrix W:\n', w)  # 射影行列w(4行×2列)を表示

# データセットを射影行列wで変換
X_train_pca = X_train_std.dot(w)

# 変換したデータをプロット
colors = ['r', 'b', 'g']
markers = ['s', 'x', 'o']

for l, c, m in zip(np.unique(y_train), colors, markers):
    plt.scatter(X_train_pca[y_train == l, 0], 
                X_train_pca[y_train == l, 1], 
                c=c, label=l, marker=m)
plt.xlabel('PC 1')
plt.ylabel('PC 2')
plt.legend(loc='lower left')
plt.tight_layout()
plt.show()
```
![スクリーンショット 2019-09-29 19.59.00.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/209705/7476deed-1aec-9ee4-8bf5-51f6b091605a.png)
上位2つの固有ベクトルを抽出し、射影行列w(4行×2列）を作成し、トレイニングデータと射影行列の内積を取ると、4次元のトレイニングデータが2次元に変換されます。

2次元に変換されたデータをみると、クラス0は明確に線形分離できそうですが、クラス1と2は少し混じってしまっている感じですね。

#4.主成分分析（sklearnライブラリー）
さて、sklearnには、主成分分析(PCA)がライブラリーにあるので、テキストのコードにある様に、これを使ってwineデータを主成分分析し、決定領域を表示をしてみましょう。

```python
from sklearn.decomposition import PCA
from sklearn.linear_model import LogisticRegression
from matplotlib.colors import ListedColormap

# 決定領域表示関数
def plot_decision_regions(X, y, classifier, resolution=0.02):

    # setup marker generator and color map
    markers = ('s', 'x', 'o', '^', 'v')
    colors = ('red', 'blue', 'lightgreen', 'gray', 'cyan')
    cmap = ListedColormap(colors[:len(np.unique(y))])

    # plot the decision surface
    x1_min, x1_max = X[:, 0].min() - 1, X[:, 0].max() + 1
    x2_min, x2_max = X[:, 1].min() - 1, X[:, 1].max() + 1
    xx1, xx2 = np.meshgrid(np.arange(x1_min, x1_max, resolution),
                           np.arange(x2_min, x2_max, resolution))
    Z = classifier.predict(np.array([xx1.ravel(), xx2.ravel()]).T)
    Z = Z.reshape(xx1.shape)
    plt.contourf(xx1, xx2, Z, alpha=0.4, cmap=cmap)
    plt.xlim(xx1.min(), xx1.max())
    plt.ylim(xx2.min(), xx2.max())

    # plot class samples
    for idx, cl in enumerate(np.unique(y)):
        plt.scatter(x=X[y == cl, 0], 
                    y=X[y == cl, 1],
                    alpha=0.6, 
                    c=cmap(idx),
                    edgecolor='black',
                    marker=markers[idx], 
                    label=cl)

# sklearnによるPCA      
pca = PCA(n_components=2)
X_train_pca = pca.fit_transform(X_train_std)
X_test_pca = pca.transform(X_test_std)
lr = LogisticRegression(solver='lbfgs', multi_class='auto')  # 引数追加
lr = lr.fit(X_train_pca, y_train)

# 学習データの決定領域表示
plot_decision_regions(X_train_pca, y_train, classifier=lr)
plt.title('PCA for traning_data')
plt.xlabel('PC 1')
plt.ylabel('PC 2')
plt.legend(loc='lower left')
plt.tight_layout()
plt.show()

# テストデータの決定領域表示
plot_decision_regions(X_test_pca, y_test, classifier=lr)
plt.title('PCA for test_data')
plt.xlabel('PC 1')
plt.ylabel('PC 2')
plt.legend(loc='lower left')
plt.tight_layout()
plt.show()

# 主成分の寄与率表示
pca = PCA(n_components=None)
X_train_pca = pca.fit_transform(X_train_std)
print('\npca.explained_variance_ratio_\n%s' %pca.explained_variance_ratio_)
```
![スクリーンショット 2019-09-29 23.55.57.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/209705/08c83bc7-2b6f-9295-6074-c496c28ae8ab.png)
sklearn のライブラリーを使うと簡単ですね。なお、sklearnのバージョン（私が使っているのは 0.21.3）の関係で、テキストのコードのままだとワーニングが出るので、lr = LogisticRegression()に引数を追加しています。