Python 機械学習プログラミングの備忘録その７（カーネル主成分分析）

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

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

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

#1.半月形データ
どう見ても線形分離は困難なデータセットの例として、半月データが登場します。これを主成分分析（PCA）で無理やり分離を試みます。

```python
import numpy as np
import matplotlib.pyplot as plt
from sklearn.decomposition import PCA
from scipy.spatial.distance import pdist, squareform
from scipy import exp
from scipy.linalg import eigh
from sklearn.datasets import make_moons
from sklearn.datasets import make_circles

# 2つの半月形データのプロット    
X, y = make_moons(n_samples=100, random_state=123)
plt.scatter(X[y == 0, 0], X[y == 0, 1], color='red', marker='^', alpha=0.5)
plt.scatter(X[y == 1, 0], X[y == 1, 1], color='blue', marker='o', alpha=0.5)
plt.title('moons_data')
plt.tight_layout()
plt.show()

# PCAによる分離
scikit_pca = PCA(n_components=2)
X_spca = scikit_pca.fit_transform(X)
fig, ax = plt.subplots(nrows=1, ncols=2, figsize=(7, 3))
ax[0].scatter(X_spca[y == 0, 0], X_spca[y == 0, 1],
              color='red', marker='^', alpha=0.5)
ax[0].scatter(X_spca[y == 1, 0], X_spca[y == 1, 1],
              color='blue', marker='o', alpha=0.5)
ax[1].scatter(X_spca[y == 0, 0], np.zeros((50, 1)) + 0.02,
              color='red', marker='^', alpha=0.5)
ax[1].scatter(X_spca[y == 1, 0], np.zeros((50, 1)) - 0.02,
              color='blue', marker='o', alpha=0.5)
ax[0].set_title('PCA for monns_data')
ax[0].set_xlabel('PC1')
ax[0].set_ylabel('PC2')
ax[1].set_title('PCA for moons_data')
ax[1].set_ylim([-1, 1])
ax[1].set_yticks([])
ax[1].set_xlabel('PC1')
plt.tight_layout()
plt.show()
```
![スクリーンショット 2019-09-30 15.36.43.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/209705/eded5e5c-6574-ac41-1c08-1ad17701ca50.png)
半月データの形が凄いですね。もうどうやっても、線形分離させないぞという形をしています。当然のことながら、主成分分析（PCA）では全く歯が立ちません。

続いて、カーネル主成分分析です。こちらは、gammaの値によって、分離具合が異なります。テキストでは、gamma=15になっていますが、1〜15まで2刻みで動かして、分離具合を連続的に見てみましょう。

```python
# カーネPCA関数
def rbf_kernel_pca(X, gamma, n_components):
    sq_dists = pdist(X, 'sqeuclidean')
    mat_sq_dists = squareform(sq_dists)
    K = exp(-gamma * mat_sq_dists)
    N = K.shape[0]
    one_n = np.ones((N, N)) / N
    K = K - one_n.dot(K) - K.dot(one_n) + one_n.dot(K).dot(one_n)
    eigvals, eigvecs = eigh(K)
    eigvals, eigvecs = eigvals[::-1], eigvecs[:, ::-1]
    X_pc = np.column_stack(list(eigvecs[:, i]
                            for i in range(n_components))) # list 追加
    return X_pc

# カーネルPCAによる分離
for i in range(1,16,2):
    X_kpca = rbf_kernel_pca(X, gamma=i, n_components=2)  
    fig, ax = plt.subplots(nrows=1,ncols=2, figsize=(7,3))
    ax[0].scatter(X_kpca[y==0, 0], X_kpca[y==0, 1], 
                color='red', marker='^', alpha=0.5)
    ax[0].scatter(X_kpca[y==1, 0], X_kpca[y==1, 1],
                color='blue', marker='o', alpha=0.5)
    ax[1].scatter(X_kpca[y==0, 0], np.zeros((50,1))+0.02, 
                color='red', marker='^', alpha=0.5)
    ax[1].scatter(X_kpca[y==1, 0], np.zeros((50,1))-0.02,
                color='blue', marker='o', alpha=0.5)
    ax[0].set_title('gamma'+str(i))
    ax[0].set_xlabel('PC1')
    ax[0].set_ylabel('PC2')
    ax[1].set_title('gamma'+str(i))
    ax[1].set_ylim([-1, 1])
    ax[1].set_yticks([])
    ax[1].set_xlabel('PC1')
    plt.tight_layout()
    plt.show()
```
![スクリーンショット 2019-09-30 15.42.15.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/209705/2c432fdc-4596-9e35-85f1-39865ac6c0fc.png)
![スクリーンショット 2019-09-30 15.43.07.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/209705/2cbfcd55-bc51-829d-5053-dbd3e1763dde.png)
![スクリーンショット 2019-09-30 15.43.25.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/209705/66952631-43ff-1e8d-920e-7df625537517.png)
おーっ、これは凄い！ gammaの値が変化するにつれて、赤と青の紐が解けて行く様な感じです。なんか、新海誠の映画に出て来そうな。。。

なお、numpyのバージョン（私は1.17.1)の関係で、X_pc = np.column_stack()のところでワーニングが出たので、listコマンドを追加しています。

#2.同心円データ
続いて、これはもう線形分離は絶対させない最終兵器とも言うべき同心円データです。一応、主成分分析もやっていますが。。。

```python
# 同心円データをプロット
X, y = make_circles(n_samples=1000, random_state=123, noise=0.1, factor=0.2)
plt.scatter(X[y == 0, 0], X[y == 0, 1], color='red', marker='^', alpha=0.5)
plt.scatter(X[y == 1, 0], X[y == 1, 1], color='blue', marker='o', alpha=0.5)
plt.title('circles_data')
plt.tight_layout()
plt.show()

# PCAによる分離
scikit_pca = PCA(n_components=2)
X_spca = scikit_pca.fit_transform(X)
fig, ax = plt.subplots(nrows=1, ncols=2, figsize=(7, 3))
ax[0].scatter(X_spca[y == 0, 0], X_spca[y == 0, 1],
              color='red', marker='^', alpha=0.5)
ax[0].scatter(X_spca[y == 1, 0], X_spca[y == 1, 1],
              color='blue', marker='o', alpha=0.5)
ax[1].scatter(X_spca[y == 0, 0], np.zeros((500, 1)) + 0.02,
              color='red', marker='^', alpha=0.5)
ax[1].scatter(X_spca[y == 1, 0], np.zeros((500, 1)) - 0.02,
              color='blue', marker='o', alpha=0.5)
ax[0].set_title('PCA for circles_data')
ax[0].set_xlabel('PC1')
ax[0].set_ylabel('PC2')
ax[1].set_title('PCA for circles_data')
ax[1].set_ylim([-1, 1])
ax[1].set_yticks([])
ax[1].set_xlabel('PC1')
plt.tight_layout()
plt.show()
```
![スクリーンショット 2019-09-30 15.52.09.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/209705/4bd41d14-8ef9-d1e1-46ba-827d4ce992b1.png)
これはどう考えても無理でしょ！主成分分析(PCA)では手も足も出ない、お手上げ状態です。

さて、カーネル主成分分析は、これを一体どう料理するのか。また、gammaを1から15まで、2刻みで変化させてみます。

```python
# カーネルPCAによる分離
for i in range(1, 16, 2): 
    X_kpca = rbf_kernel_pca(X, gamma=i, n_components=2)
    fig, ax = plt.subplots(nrows=1, ncols=2, figsize=(7, 3))
    ax[0].scatter(X_kpca[y == 0, 0], X_kpca[y == 0, 1],
                  color='red', marker='^', alpha=0.5)
    ax[0].scatter(X_kpca[y == 1, 0], X_kpca[y == 1, 1],
                  color='blue', marker='o', alpha=0.5)
    ax[1].scatter(X_kpca[y == 0, 0], np.zeros((500, 1)) + 0.02,
                  color='red', marker='^', alpha=0.5)
    ax[1].scatter(X_kpca[y == 1, 0], np.zeros((500, 1)) - 0.02,
                  color='blue', marker='o', alpha=0.5)
    ax[0].set_title('gamma'+str(i))
    ax[0].set_xlabel('PC1')
    ax[0].set_ylabel('PC2')
    ax[1].set_title('gamma'+str(i))
    ax[1].set_ylim([-1, 1])
    ax[1].set_yticks([])
    ax[1].set_xlabel('PC1')
    plt.tight_layout()
    plt.show()
```
![スクリーンショット 2019-09-30 15.56.13.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/209705/c4332118-4440-9844-2a19-3edf2fdae3b0.png)
![スクリーンショット 2019-09-30 15.56.33.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/209705/07d04bb2-341a-64f3-7df4-6acc622c57c3.png)
![スクリーンショット 2019-09-30 15.56.51.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/209705/ee6625f8-a57b-9e27-c086-06cfbb9e2fd5.png)
これはもの凄い！なんと細胞分裂しちゃいました。うむむ、これなら線形分離可能です。まさしくなんでも有りだぞ、カーネル主成分分析！

#4.新しいデータ点の射影
カーネルPCA主成分分析で、トレーニングに使った以外のデータを分析する方法がテキストに出ています。実行すると、

```python
# カーネル行列関数（固有値を返す様に修正）
def rbf_kernel_pca(X, gamma, n_components):
    sq_dists = pdist(X, 'sqeuclidean')
    mat_sq_dists = squareform(sq_dists)
    K = exp(-gamma * mat_sq_dists)
    N = K.shape[0]
    one_n = np.ones((N, N)) / N
    K = K - one_n.dot(K) - K.dot(one_n) + one_n.dot(K).dot(one_n)
    eigvals, eigvecs = eigh(K)
    eigvals, eigvecs = eigvals[::-1], eigvecs[:, ::-1]
    alphas = np.column_stack(list(eigvecs[:, i]
                              for i in range(n_components)))  # list 追加
    lambdas = [eigvals[i] for i in range(n_components)]
    return alphas, lambdas


X, y = make_moons(n_samples=100, random_state=123)
alphas, lambdas = rbf_kernel_pca(X, gamma=15, n_components=1)

# 26番目の点
x_new = X[25]
print('x_new = ', x_new)
x_proj = alphas[25] 
print('x_pro = ', x_proj)

# 射影関数
def project_x(x_new, X, gamma, alphas, lambdas):
    pair_dist = np.array([np.sum((x_new - row)**2) for row in X])
    k = np.exp(-gamma * pair_dist)
    return k.dot(alphas / lambdas)

# 新たなデータの射影
x_reproj = project_x(x_new, X, gamma=15, alphas=alphas, lambdas=lambdas)
x_reproj

# 射影した点をプロット 
plt.scatter(alphas[y == 0, 0], np.zeros((50)),
            color='red', marker='^', alpha=0.5)
plt.scatter(alphas[y == 1, 0], np.zeros((50)),
            color='blue', marker='o', alpha=0.5)
plt.scatter(x_proj, 0, color='black',
            label='original projection of point X[25]', marker='^', s=100)
plt.scatter(x_reproj, 0, color='green',
            label='remapped point X[25]', marker='x', s=500)
plt.legend(scatterpoints=1)
plt.tight_layout()
plt.show()
```
![スクリーンショット 2019-09-30 16.22.24.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/209705/a465d79e-5cd8-5f15-67db-6084729c0a28.png)
正しく射影されましたね。

さて、最後に、ssklearnのライブラリーにあるカーネル主成分分析のコードを動かしましょう。

```python
# sklearn ライブラリー　KernelPCA
from sklearn.decomposition import KernelPCA

X, y = make_moons(n_samples=100, random_state=123)
scikit_kpca = KernelPCA(n_components=2, kernel='rbf', gamma=15)
X_skernpca = scikit_kpca.fit_transform(X)

plt.scatter(X_skernpca[y == 0, 0], X_skernpca[y == 0, 1],
            color='red', marker='^', alpha=0.5)
plt.scatter(X_skernpca[y == 1, 0], X_skernpca[y == 1, 1],
            color='blue', marker='o', alpha=0.5)

plt.xlabel('PC1')
plt.ylabel('PC2')
plt.tight_layout()
plt.show()
```
![スクリーンショット 2019-09-30 16.24.53.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/209705/80edd7ba-99d9-45d6-ecb9-4295191ff336.png)
ライブラリーを使うと簡単！