次元削減でデータを圧縮する
データ圧縮は機械学習の重要なテーマの一つであり、データ量の多い現代のテクノロジにおいてデータ格納と分析に役立つ。
4章では次元削減のうちの特徴選択を行なった
5章では同じく次元削減の特徴抽出を行う
特徴選択と同様特量抽出を使用すればデータセット内の特徴量の個数を減らすことができる
特徴選択では、元の特徴量は変換されることなくそのままの形で維持されていたが、
特徴抽出ではデータが新しい特徴空間に変換、または射影される
本章では次元抽出の中でも以下の三つを取り上げる
①主成分分析（PCA)
②線形判別分析（LDA）
③カーネル主成分分析
①主成分分析による教師なし次元削減
pacは様々な分野にわたって広く使用されている教師なし線形変換法である
【活用例】
・探索的データ解析
・株取引での信号のノイズ除去
・ゲノムデータや遺伝子発現量の解析
特徴量空間において，意味のない（相関の高い）特徴量の組合わせを削除する．
pcaは二つの方法がある
・分散最大基準←今回はこっち
・平均二乗誤差最小基準
流れとしては
1,分散が大きい(違いがよく表されている)軸を見つける
2,その軸に直行する分散が大きい軸を見つける
3,2をk回繰り返す
元のd次元のデータを新しいk次元の部分空間に変換すると、最初の主成分の分散は最大となる
結果として生じるすべての主成分の分散が最大となるのは、他の成分と相関がない場合となる
※PCAの方向はデータのスケーリングに対して非常に敏感である。そのため、最初に特徴量を標準化する必要がある
http://www.me.cs.scitec.kobe-u.ac.jp/~takigu/jugyou/shingou/0210/ref/r8_pca.pdf
http://www.slideshare.net/takanoriogata1121/10pca-49324044?next_slideshow=1
実際に行うステップは以下の様になる
1,d次元のデータセットを標準化する
2,標準化したデータセットの共分散行列を作成する
3,共分散行列を固有ベクトルと固有値に分解する
4,最も大きいk個の固有値に対応するk個の固有ベクトルを選択する
5,上位k個の固有ベクトルから射影行列Wを作成する
6,射影行列Wを使ってd次元の入力データセットXを変換し、新しい次元の特徴部分空間を取得する
まずWINEデータセットをロードする
import pandas as pd
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']
df_wine.head()
d次元のデータセットを標準化する
トレーングデータとテストデータに分割し、分散が1となるよう標準化する
if Version(sklearn_version) < '0.18':
from sklearn.cross_validation import train_test_split
else:
from sklearn.model_selection import train_test_split
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, random_state=0)
from sklearn.preprocessing import StandardScaler
sc = StandardScaler()
X_train_std = sc.fit_transform(X_train)
X_test_std = sc.transform(X_test)
2,標準化したデータセットの共分散行列を作成する
二つの特徴量xjとxkの間の共分散は以下の式を使って計算できる
σ_{jk} = \frac{1}{n}\sum_{i=1}^{n}(x_j^i-μ_j)(x_k^i-μ_k) \tag{5.1.2}
Σ=
\begin{pmatrix}
σ^2_1 & σ_12 & σ_13\\
σ_21 & σ^2_2 & σ_23\\
σ_31 & σ_32 & σ^2_3
\end{pmatrix}
\tag{5.1.3}
Σ_v = λ_v \tag{5.1.4}
こんな計算を手で計算するのは大変なのでNumPyのlinalg.eig関数を使ってWineデータセットの共分散行列の固有対を取得する
import numpy as np
cov_mat = np.cov(X_train_std.T)
eigen_vals, eigen_vecs = np.linalg.eig(cov_mat)
print('\nEigenvalues \n%s' % eigen_vals)
3,共分散行列を固有ベクトルと固有値に分解する
データセットを新しい特徴部分空間に圧縮するという方法でデータセットの次元を削減したい
そのため、データに含まれる大半の情報を含んでいる固有ベクトルだけを選択する
固有値は固有ベクトルの大きさを表すため、大きいものから順に降順で並び変える必要がある
ここでは、固有値の分散説明率をプロットしてみる
分散説明率とは、固有値の合計に対する固有値の割合のことである
\frac{λ_j}{Σ^d_{j=1}λ_j}
\tag{5.1.5}
分散説明率の累積和はNumpyのcumsum関数を使って計算できる
計算した値のプロットはmatplotlibのstep関数を使用する
tot = sum(eigen_vals)
var_exp = [(i / tot) for i in sorted(eigen_vals, reverse=True)]
cum_var_exp = np.cumsum(var_exp)
import matplotlib.pyplot as plt
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 components')
plt.legend(loc='best')
plt.tight_layout()
# plt.savefig('./figures/pca1.png', dpi=300)
plt.show()
一つ目と二つ目の主成分だけで全体の60%ほどとなることがわかる
4,最も大きいk個の固有値に対応するk個の固有ベクトルを選択する
固有値の大きいものから順に固有対を並べ替え、選択された固有ベクトルから射影行列を生成する
まず、固有値の大きいものから順に固有対を並び変える
# Make a list of (eigenvalue, eigenvector) tuples
eigen_pairs = [(np.abs(eigen_vals[i]), eigen_vecs[:, i])
for i in range(len(eigen_vals))]
# Sort the (eigenvalue, eigenvector) tuples from high to low
eigen_pairs.sort(key=lambda k: k[0], reverse=True)
# Note: I added the `key=lambda k: k[0]` in the sort call above
# just like I used it further below in the LDA section.
# This is to avoid problems if there are ties in the eigenvalue
# arrays (i.e., the sorting algorithm will only regard the
# first element of the tuples, now).
次に、最も二つの固有値に対する二つの固有ベクトルを集める
これによって分散の60%をおさえられる
※実務では、計算効率と分類器の性能のバランスを見ながら、主成分の個数を決定する必要がある
5,上位k個の固有ベクトルから射影行列Wを作成する
二つの固有ベクトルから13×2次元の射影行列Wが生成される
w = np.hstack((eigen_pairs[0][1][:, np.newaxis],
eigen_pairs[1][1][:, np.newaxis]))
print('Matrix W:\n', w)
6,射影行列Wを使ってd次元の入力データセットXを変換し、新しい次元の特徴部分空間を取得する
行列の積を計算することで、124×13次元のトレーニングデータセット全体を2つの主成分に変換できる
X' = XW \tag{5.1.7}
X_train_pca = X_train_std.dot(w)
最後に、124×2次元の行列として格納された変換後のwineトレーニングセットを、二次元の散布図としてプロットしてみる
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.savefig('./figures/pca2.png', dpi=300)
plt.show()
scikitlearnを使ってみる
1,2,3をまとめてやってみると
from sklearn.decomposition import PCA
pca = PCA()
X_train_pca = pca.fit_transform(X_train_std)
pca.explained_variance_ratio_
4,5,6をまとめてやると
pca = PCA(n_components=2)
X_train_pca = pca.fit_transform(X_train_std)
X_test_pca = pca.transform(X_test_std)
plt.scatter(X_train_pca[:, 0], X_train_pca[:, 1])
plt.xlabel('PC 1')
plt.ylabel('PC 2')
plt.show()
二つの主成分軸に削減されたトレーニングモデルの決定領域をみてみる
from sklearn.linear_model import LogisticRegression
lr = LogisticRegression()
lr = lr.fit(X_train_pca, y_train)
plot_decision_regions(X_train_pca, y_train, classifier=lr)
plt.xlabel('PC 1')
plt.ylabel('PC 2')
plt.legend(loc='lower left')
plt.tight_layout()
# plt.savefig('./figures/pca3.png', dpi=300)
plt.show()
線形判別分類による教師ありデータ圧縮
線形判別分析は、特徴抽出手法の一つである。
LDAの基本的な考え方はPCAとよく似ているが、PCAはデータセットにおいて分散がもっとも大きい直行成分軸を見つけ出そうとするのに対して、LDAはクラスの分離を最適化する特徴部分空間を見つけ出そうとする。
LDAとPCAはどちらも線形変換方であり、data setの次元数を減らすために使用できる。
LDAとPCAはどちらも線形変換法であり、データセットの次元数を減らすために使用できる。
PCAが教師なしのアルゴリズムであるのに対し、LDAは教師ありのアルゴリズムである。
このため、分類タスクの特徴抽出の手段としては、LDAの方が優れている。
LDAではデータが正規分布に従っていることが前提となる。また、クラスの共分散行列がまったくおなじであることと、特徴量が統計的にみて互いに独立していることも前提となる。
※ただし、前提条件を多少満たしていなくても、次元削減の手段としてのLDAはそれなりにうまくくようです
LDAの手順
1,d次元のデータセットを標準化する
2,クラス毎にd次元の平均ベクトルを計算する
3,クラス間変動行列Sbと、クラス内変動係数Swを生成する
4,行列Sx^-1Sbの固有ベクトルと対応する固有値を計算する
5,d×k次元の変換行列Wを生成するために、最も大きいk個の固有値に対応するk個の固有ベクトルを選択する。固有ベクトルは、この行列の列である。
6,変換行列Wを使ってサンプルを新しい特徴部分空間へ射影する
1,d次元のデータセットを標準化する
wineデータセットの特徴量は、PCAに関する最初の節ですでに標準化している。このため、最初のステップを飛ばして、平均ベクトルの計算に進むことができる。
2,クラス毎にd次元の平均ベクトルを計算する
各平均ベクトルmiんは、クラスiのサンプルに関する平均特徴量の値μmが格納される
np.set_printoptions(precision=4)
mean_vecs = []
for label in range(1, 4):
mean_vecs.append(np.mean(X_train_std[y_train == label], axis=0))
print('MV %s: %s\n' % (label, mean_vecs[label - 1]))
3,クラス間変動行列Sbと、クラス内変動係数Swを生成する
平均ベクトルを使ってクラス内変動行列Swを計算する方法は以下のようになる
【式5.2.3】
式5.2.3の右辺は、個々のクラスiについて変動行列Sjを合計することによって計算される
【式5.2.4】
d = 13 # number of features
S_W = np.zeros((d, d))
for label, mv in zip(range(1, 4), mean_vecs):
class_scatter = np.zeros((d, d)) # scatter matrix for each class
for row in X_train_std[y_train == label]:
row, mv = row.reshape(d, 1), mv.reshape(d, 1) # make column vectors
class_scatter += (row - mv).dot((row - mv).T)
S_W += class_scatter # sum class scatter matrices
print('Within-class scatter matrix: %sx%s' % (S_W.shape[0], S_W.shape[1]))
変動行列を計算するときには、トレーニングデータセットにおいてクラスラベルが一様に分布していることが前提となる。
そこで調べてみる
print('Class label distribution: %s'
% np.bincount(y_train)[1:])
Class label distribution: [40 49 35]
しかし、クラスラベル個数を出力すると前提を満たしていないことがわかる
従って、スケーリングが必要となる
（変動行列をクラスのサンプルの個数で割る）
【式5.2.5】
d = 13 # number of features
S_W = np.zeros((d, d))
for label, mv in zip(range(1, 4), mean_vecs):
class_scatter = np.cov(X_train_std[y_train == label].T)
S_W += class_scatter
print('Scaled within-class scatter matrix: %sx%s' % (S_W.shape[0],
S_W.shape[1]))
スケーリングされたクラス内変動行列を計算した後は、
クラス間変動行列Sbの計算に進む
【式5.2.6】
mean_overall = np.mean(X_train_std, axis=0)
d = 13 # number of features
S_B = np.zeros((d, d))
for i, mean_vec in enumerate(mean_vecs):
n = X_train[y_train == i + 1, :].shape[0]
mean_vec = mean_vec.reshape(d, 1) # make column vector
mean_overall = mean_overall.reshape(d, 1) # make column vector
S_B += n * (mean_vec - mean_overall).dot((mean_vec - mean_overall).T)
print('Between-class scatter matrix: %sx%s' % (S_B.shape[0], S_B.shape[1]))
4,行列Sx^-1Sbの固有ベクトルと対応する固有値を計算する
行列Sw~-1Sbの一般化された固有値問題を解く
eigen_vals, eigen_vecs = np.linalg.eig(np.linalg.inv(S_W).dot(S_B))
固有値を計算した後は、固有値を大きいものから降順で並べ替えることができる
# Make a list of (eigenvalue, eigenvector) tuples
eigen_pairs = [(np.abs(eigen_vals[i]), eigen_vecs[:, i])
for i in range(len(eigen_vals))]
# Sort the (eigenvalue, eigenvector) tuples from high to low
eigen_pairs = sorted(eigen_pairs, key=lambda k: k[0], reverse=True)
# Visually confirm that the list is correctly sorted by decreasing eigenvalues
print('Eigenvalues in decreasing order:\n')
for eigen_val in eigen_pairs:
print(eigen_val[0])
クラス間変動行列Sbは、階数1以下のc個の行列を合計したものである。
5,d×k次元の変換行列Wを生成するために、最も大きいk個の固有値に対応するk個の固有ベクトルを選択する。固有ベクトルは、この行列の列である。
線形判別によってクラスの判別情報がどの程度細くされるのかを測定するため、PCAの節で作成した分散説明のグラフと同様に、固有値を減らしながら線形判別をプロットしてみる
※ここではクラスの判別情報を「判別性」と呼ぶ
tot = sum(eigen_vals.real)
discr = [(i / tot) for i in sorted(eigen_vals.real, reverse=True)]
cum_discr = np.cumsum(discr)
plt.bar(range(1, 14), discr, alpha=0.5, align='center',
label='individual "discriminability"')
plt.step(range(1, 14), cum_discr, where='mid',
label='cumulative "discriminability"')
plt.ylabel('"discriminability" ratio')
plt.xlabel('Linear Discriminants')
plt.ylim([-0.1, 1.1])
plt.legend(loc='best')
plt.tight_layout()
# plt.savefig('./figures/lda1.png', dpi=300)
plt.show()
結果のグラフからわかるように最初の二つの線形判別はwineトレーニングデータセット内の情報をほぼ100%補足している
なので二つの固有ベクトルを列方向に並べて変換行列Wを作成してみる
w = np.hstack((eigen_pairs[0][1][:, np.newaxis].real,
eigen_pairs[1][1][:, np.newaxis].real))
print('Matrix W:\n', w)
6,変換行列Wを使ってサンプルを新しい特徴部分空間へ射影する
作成した変換行列Wに行列を掛けることにより、トレーニングデータセットを変換してみる
X' = XW
X_train_lda = 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_lda[y_train == l, 0] * (-1),
X_train_lda[y_train == l, 1] * (-1),
c=c, label=l, marker=m)
plt.xlabel('LD 1')
plt.ylabel('LD 2')
plt.legend(loc='lower right')
plt.tight_layout()
# plt.savefig('./figures/lda2.png', dpi=300)
plt.show()
結果として得られたグラフをみると、三つのクラスが完全に線形分離となる
scikit-learnによるLDA
if Version(sklearn_version) < '0.18':
from sklearn.lda import LDA
else:
from sklearn.discriminant_analysis import LinearDiscriminantAnalysis as LDA
lda = LDA(n_components=2)
X_train_lda = lda.fit_transform(X_train_std, y_train)
以上でこの節でやったことは終了
ではロジスティック回帰の分類器がLDAによって変換された低次元のトレーニングデータセットをどのように処理するのかみてみる
from sklearn.linear_model import LogisticRegression
lr = LogisticRegression()
lr = lr.fit(X_train_lda, y_train)
plot_decision_regions(X_train_lda, y_train, classifier=lr)
plt.xlabel('LD 1')
plt.ylabel('LD 2')
plt.legend(loc='lower left')
plt.tight_layout()
# plt.savefig('./images/lda3.png', dpi=300)
plt.show()
【画像】
グラフをみるとロジスティック回帰モデルがクラス2のサンプルの一つを誤分類していることがわかる
次にテストデータの結果を見てみる
X_test_lda = lda.transform(X_test_std)
plot_decision_regions(X_test_lda, y_test, classifier=lr)
plt.xlabel('LD 1')
plt.ylabel('LD 2')
plt.legend(loc='lower left')
plt.tight_layout()
# plt.savefig('./images/lda4.png', dpi=300)
plt.show()
【画像】
ロジスティック回帰では、wineのデータセットの元の13個の特徴量ではなく、二次元の特徴部分空間だけを使用することによりテストデータセットのサンプルを分類できることがわかる
カーネル主成分分析を使った非線形写像
この節では、カーネル化されたPCA、カーネルPCAを使用して、線形に分離できないデータを変換し、線形分類器に適した新しい低次元の部分空間へ射影する方法を学ぶ
すみませんやってる事の意味がわかりませんでした
例1
こんな形のデータセット
これは線形分離はできない
ここではカーネルPCAを用いて半月を「展開」し、データセットを線形分類器に適した入力にする
しかし、その前に標準のPCAを使ってやってみる
標準のPCAを使って変換したものでは線形分類器が十分な性能が出せないことがわかる
では次にカーネルPCAを使ってみると
線形分類できるようになった
次に同心円の分離を行う
これも線形分離はできない
だが一応PCAをあててみる
やはりうまくいかない
そこでカーネルPCAを行なってみる
できた
scikit-learnでやってみる
最初にやった半月型の例
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.savefig('./figures/scikit_kpca.png', dpi=300)
plt.show()
先ほどの結果と一致している。
以上です。