機械学習
1. 線形回帰モデル
1.1 要点
回帰問題とは、ある入力から出力を予測する問題のことである。直線で予測する場合を線形回帰、曲線で予測する場合は非線形回帰と呼ぶ。
線形回帰モデルは回帰問題を解くための教師あり学習の1つであり、
入力とパラメータの線形結合を出力するモデル
入力$\boldsymbol{x}$、パラメータ$\boldsymbol{w}$を、
\boldsymbol{x}=(x_1,x_2,_\cdots,x_m)^T\\
\boldsymbol{w}=(w_1,w_2,_\cdots,w_m)^T
としたとき、出力$\hat{y}$(予測値)は、
\hat{y}=\boldsymbol{w}^T\boldsymbol{x}+w_0=\sum_{j=1}^mw_jx_j+w_0
と表される。
説明変数$x_i$が1つのとき単回帰と呼び、2つ以上のとき重回帰と呼ぶ。
次のような連立方程式を考える。
y_1=w_0+w_1x_{11}+w_2x_{12}+\cdots+w_mx_{1m}+\epsilon_1\\
y_2=w_0+w_1x_{21}+w_2x_{22}+\cdots+w_mx_{2m}+\epsilon_2\\
\vdots\\
y_n=w_0+w_1x_{n1}+w_2x_{n2}+\cdots+w_mx_{nm}+\epsilon_n
これを行列表現で表すと以下のようになる。
\boldsymbol{y}=\boldsymbol{X}\boldsymbol{w}+\boldsymbol{\epsilon}\\
\boldsymbol{X}=(\boldsymbol{x_1},\cdots,\boldsymbol{x_n})^T\\
\boldsymbol{y}=(y_1,\cdots,y_n)^T\\
\boldsymbol{x_i}=(1,x_{i1},\cdots,x_{im})^T\\
\boldsymbol{\epsilon}=(\epsilon_1,\cdots,\epsilon_n)^T\\
\boldsymbol{w}=(w_0,w_1,\cdots,w_m)^T
このとき$\boldsymbol{X}$は計画行列と呼ぶ。
最小二乗法
パラメータ$\boldsymbol{w}$は最小二乗法で決められる。
**平均二乗誤差（残差平方和）**とは、データとモデル出力の二乗誤差の和であり、小さいほど直線とデータの距離が近い。
MSE_{train}=\frac{1}{n_{train}}\sum_{i=1}^{n_{train}}(\hat{y}_i^{(train)}-y_i^{(train)})^2
ここで、最小二乗法による$\boldsymbol{w}$の推定を行う。
MSEを最小化するには、$\boldsymbol{w}$に関して微分したものが0になる点を求めればよい。
\frac{\partial{}}{\partial \boldsymbol{w}}\biggl\{\frac{1}{n_{train}}\sum_{i=1}^{n_{train}}(\hat{y}_i^{(train)}-y_i^{(train)})^2 \biggr\}=0
\frac{1}{n_{train}}\frac{\partial{}}{\partial \boldsymbol{w}}\bigl\{(\boldsymbol{X}^{(train)}\boldsymbol{w}-\boldsymbol{y}^{(train)})^T(\boldsymbol{X}^{(train)}\boldsymbol{w}-\boldsymbol{y}^{(train)}) \bigr\}=0
2\boldsymbol{X}^{(train)T}\boldsymbol{X}^{(train)}\boldsymbol{w}-2\boldsymbol{X}^{(train)T}\boldsymbol{y}^{(train)}=0
\hat{\boldsymbol{w}}=(\boldsymbol{X}^{(train)T}\boldsymbol{X}^{(train)})^{-1}\boldsymbol{X}^{(train)T}\boldsymbol{y}^{(train)}
予測値$\hat{\boldsymbol{y}}$は
\begin{align}
\hat{\boldsymbol{y}}&=\boldsymbol{X}\hat{\boldsymbol{w}}\\
&=\boldsymbol{X}・(\boldsymbol{X}^{(train)T}\boldsymbol{X}^{(train)})^{-1}\boldsymbol{X}^{(train)T}\boldsymbol{y}^{(train)}
\end{align}
となる。
1.2 実装
shikit-learnを用いて線形回帰モデルの実装を行う。
ここで使用するデータはボストン近郊の住宅価格のデータセットであり、これはshikit-learnで用意されているものである。
from sklearn.datasets import load_boston
import pandas as pd
import numpy as np
boston = load_boston()
df = pd.DataFrame(boston.data, columns=boston.feature_names)
df['PRICE'] = boston.target
df.head()
ここでは、目的変数は住宅価格であり、説明変数として、犯罪発生率(CRIM)、一戸当たりの平均部屋数(RM)を使用する。
data = df.loc[:, ['CRIM', 'RM']].values
target = df.loc[:,['PRICE']].values
モデルの作成をする。
線形回帰モデルはshikit-learnのLinearRegressionを使って実装できる。
from sklearn.linear_model import LinearRegression
model = LinearRegression()
model.fit(data, target)
そして「predict」で与えた説明変数に対する予測結果を出力することができる。
model.predict([[0.2, 7]])
切片と回帰係数の確認を行う。
print('回帰係数: ', model.coef_)
print('切片: ', model.intercept_)
2. 非線形回帰モデル
2.1 要約
非線形回帰とは、
y = w_0+w_1x+w_2x^2+w_3x^3\\
y = w_0+w_1\sin{x}+w_2\cos{x}+w_3\log{x}
などのように、
線形回帰での$x$の代わりに$x$の関数$\phi (x)$を用いる回帰モデルのこと。
ただし、パラメータ$\boldsymbol{w}$については線形の関係であることは重回帰モデルと変わりはない。
基底展開法
回帰関数として、基底関数と呼ばれる既知の非線形関数とパラメータベクトルの線形結合を使用
y_i=w_0+\sum_{j=1}^mw_j\phi(\boldsymbol{x_i})+\epsilon_j
よく使われる基底関数には多項式関数やガウス基底関数などがある。
基底展開法も線形回帰と同じ枠組みで推定が可能である。
説明変数を
\boldsymbol{x}_i=(x_{i1},x_{i2},\cdots,x_{im})
非線形関数ベクトルを
\boldsymbol{\phi}(\boldsymbol{x}_i)=(\phi_1(\boldsymbol{x}_i),\phi_2(\boldsymbol{x}_i),\cdots,\phi_k(\boldsymbol{x}_i))
としたとき、計画行列は
\Phi^{(train)}=(\boldsymbol{\phi}(\boldsymbol{x}_1),\boldsymbol{\phi}(\boldsymbol{x}_2),\cdots,\boldsymbol{\phi}(\boldsymbol{x}_n))
となることから、予測値は最尤法により
\hat{\boldsymbol{y}}=\boldsymbol{\Phi}(\Phi^{(train)T}\Phi^{(train)})^{-1}\Phi^{(train)T}\boldsymbol{y}^{(train)}
と計算できる。
未学習と過学習
学習データに対して十分小さな誤差が得られないモデルを未学習と呼ぶ。
小さな誤差であるが、テストデータに対する誤差との差が大きいモデルを過学習と呼ぶ。
対策としては、
-
学習データの数を増やす
-
不要な基底関数を削除(特徴量選択)
モデルの複雑さが変化 -
正則化法を利用
モデルの複雑さに従って、値が大きくなる**正則化項(罰則項)**を課した関数を最小化
などがある。
正則化項
- L2ノルム Ridge推定量と呼ばれパラメータが0に近づくように推定される
- L1ノルム Lasso推定量と呼ばれいくつかのパラメータを正確に0に推定される
正則化パラメータは大きいほど制約面が小さくなり、パラメータが0に近くなるように推定される。
モデルが未学習か過学習しているかの判断としては、
- 訓練誤差もテスト誤差もどちらも小さい 汎化している
- 訓練誤差は小さいがテスト誤差が大きい 過学習
- 訓練誤差もテスト誤差もどちらも小さくならない 未学習
精度評価
ホールドアウト法
データを学習用データとテスト用に分割し、予測精度や誤り率を推定するために使用
手元にデータが大量にある場合を除いて、よい性能評価を与えないという欠点がある
クロスバリデーション
データを複数のグループに分割し、
それぞれのグループを検証用データ、その他のデータを学習用データとしたモデルを作成し精度評価を行う。
これらの精度の平均値をCV値と呼ぶ。
パラメータのチューニング
グリッドサーチ
- すべてのチューニングパラメータの組み合わせで評価値を算出
- 最も良い評価値を持つチューニングパラメータを持つ組み合わせを採用
2.2 実装
未学習と過学習の様子を示す。
次数が低いと表現力は低く、次数が高くなると表現力は高い分複雑なモデルとなる。
import matplotlib.pyplot as plt
import japanize_matplotlib
x = np.arange(0,10, 0.2)
def poly(x):
return 1 + 0.8*x + 0.3*x**2 - 0.04*x**3 + 3*np.random.rand(len(x))
X1 = np.array([x]).T
X2 = np.array([x, x**2, x**3]).T
X3 = np.array([x, x**2, x**3, x**4, x**5, x**6, x**7, x**8, x**9, x**10]).T
lr1 = LinearRegression()
lr2 = LinearRegression()
lr3 = LinearRegression()
lr1.fit(X1, poly(x))
lr2.fit(X2, poly(x))
lr3.fit(X3, poly(x))
plt.plot(poly(x), 'o')
plt.plot(lr1.predict(X1), color='red', label='1次')
plt.plot(lr2.predict(X2), color='green', label='3次')
plt.plot(lr3.predict(X3), color='blue', label='10次')
plt.legend();
Ridge回帰の様子を確認する。
ここでは大きな違いは見られなかった。
from sklearn.linear_model import Ridge
ridge1 = Ridge(alpha=0.1)
ridge2 = Ridge(alpha=1.0)
ridge3 = Ridge(alpha=10.0)
ridge1.fit(X3, poly(x))
ridge2.fit(X3, poly(x))
ridge3.fit(X3, poly(x))
plt.plot(poly(x), 'o')
plt.plot(ridge1.predict(X3), color='red', label='alpha=0.1')
plt.plot(ridge2.predict(X3), color='green', label='alpha=1.0')
plt.plot(ridge3.predict(X3), color='blue', label='alpha=10.0')
plt.legend();
それぞれの係数を確認する。
正則化のパラメータが高くなるにつれて係数が小さくなっていくことが確認できた。
plt.plot(ridge1.coef_, 'o', color='red', label='alpha=0.1')
plt.plot(ridge2.coef_, 'o', color='green', label='alpha=1.0')
plt.plot(ridge3.coef_, 'o', color='blue', label='alpha=10.0')
plt.plot(lr3.coef_, 'o', color='orange', label='alpha=0')
plt.legend();
plt.xlabel('xの次数');
plt.ylabel('係数');
同様にLassoでも確認を行う。
from sklearn.linear_model import Lasso
lasso1 = Lasso(alpha=0.1)
lasso2 = Lasso(alpha=1.0)
lasso3 = Lasso(alpha=10.0)
lasso1.fit(X3, poly(x))
lasso2.fit(X3, poly(x))
lasso3.fit(X3, poly(x))
plt.plot(poly(x), 'o')
plt.plot(lasso1.predict(X3), color='red', label='alpha=0.1')
plt.plot(lasso2.predict(X3), color='green', label='alpha=1.0')
plt.plot(lasso3.predict(X3), color='blue', label='alpha=10.0')
plt.legend();
基の回帰係数に比べほとんどの係数が0になっていることが確認できた。
plt.plot(lasso1.coef_, 'o', color='red', label='alpha=0.1')
plt.plot(lasso2.coef_, 'o', color='green', label='alpha=1.0')
plt.plot(lasso3.coef_, 'o', color='blue', label='alpha=10.0')
plt.plot(lr3.coef_, 'o', color='orange', label='alpha=0')
plt.legend();
plt.xlabel('xの次数');
plt.ylabel('係数');
次に、交差検証の例を示す。ボストン住宅データを使用し、5分割交差検証を行う。
from sklearn.model_selection import train_test_split
from sklearn.model_selection import KFold
from sklearn.metrics import mean_squared_error
data = df.loc[:, ['CRIM', 'RM']].values
target = df.loc[:,['PRICE']].values
x_train, x_test, y_train, y_test = train_test_split(data, target)
kf = KFold(n_splits=5)
train_loss = []
valid_loss = []
for train_index, valid_index in kf.split(x_train):
x_train_kf, x_valid_kf = x_train[train_index, :], x_train[valid_index, :]
y_train_kf, y_valid_kf = y_train[train_index, :], y_train[valid_index, :]
ridge = Ridge(alpha=0.1)
ridge.fit(x_train_kf, y_train_kf)
trainingLoss = mean_squared_error(y_train_kf, ridge.predict(x_train_kf))
train_loss.append(trainingLoss)
validLoss = mean_squared_error(y_valid_kf, ridge.predict(x_valid_kf))
valid_loss.append(validLoss)
print('train loss: ',trainingLoss, 'validation loss: ', validLoss)
最後にグリッドサーチの例を示す。
ここではRidge回帰の正則化パラメータの値のチューニングを行う。
from sklearn.model_selection import GridSearchCV
param_grid = {'alpha':[0.001, 0.01, 0.1, 1.0, 10.0, 100.0]}
# refit=Trueで訓練データ全体を使って最適な設定で自動的に再適合する
gridsearch = GridSearchCV(Ridge(), param_grid, cv=5, verbose=0, refit=True, scoring='neg_mean_squared_error')
best_model = gridsearch.fit(x_train, y_train)
print(-best_model.best_score_)
print(best_model.best_params_)
3. ロジスティック回帰モデル
3.1 要約
分類問題とは、ある入力からクラスに分類する問題のことである。
分類問題に対するアプローチ
- 識別的アプローチ
$p(C_k|\boldsymbol{x})$を直接モデル化 - 生成的アプローチ
$p(C_k)$と$p(\boldsymbol{x}|C_k)$をモデル化し、Bayesの定理を用いて$p(C_k|\boldsymbol{x})$を求める
ロジスティック回帰は識別的アプローチである。
入力$\boldsymbol{x}$、パラメータ$\boldsymbol{w}$を、
\boldsymbol{x}=(x_1,x_2,_\cdots,x_m)^T\\
\boldsymbol{w}=(w_1,w_2,_\cdots,w_m)^T
としたとき、これらの線形結合$\hat{y}$は、
\hat{y}=\boldsymbol{w}^T\boldsymbol{x}+w_0=\sum_{j=1}^mw_jx_j+w_0
となる。
この線形結合をシグモイド関数に入力することで、y=1になる確率を出力する。
シグモイド関数とは、
\sigma{(x)}=\frac{1}{1+\exp{(-ax)}}
で表される関数で、出力は必ず0-1の値となる。
シグモイド関数の微分
\sigma{(z)}=\frac{1}{1+\exp{(-z)}}
の微分は、
\frac{\partial \sigma{(z)}}{\partial z}=(-1) * (1+\exp{(-z)})^{-2}*\exp{(-z)}*(-1)\\
=\sigma{(z)}(1-\sigma{(z)})
と計算できる。
ロジステック回帰モデルは、シグモイド関数を使用して
P(Y=1|\boldsymbol{x})=\sigma(w_0+w_1x_1+\cdots+w_mx_m)
と表せる。
このときY=1になる確率が0.5以上のとき1、0.5未満のときに0とするなど判断できる。
最尤推定
y1～ynのデータが得られた際の尤度関数は
\begin{align}
P(y_1,y_2,\cdots,y_n;p)&=\prod_{i=1}^np^{y_i}(1-p)^{1-y_i}\\
&=\prod_{i=1}^n\sigma{(\boldsymbol{w^T}\boldsymbol{x_i})}^{y_i}(1-\sigma{(\boldsymbol{w^T}\boldsymbol{x_i})})^{1-y_i}\\
&=L(\boldsymbol{w})
\end{align}
となり、この尤度関数を最大化するパラメータを探索する。
対数をとることで微分の計算が可能となる。
\begin{align}
E(w_0,w_1,\cdots,w_m)&=-\log{L(\boldsymbol{w})}\\
&=\sum_{i=1}^n\{y_i\log{p_i}+(1-y_i)\log{(1-p_i)}\}
\end{align}
勾配降下法
反復学習によりパラメータを逐次的に更新する
学習率$\eta$の調整により収束のしやすさを調整
ここで対数尤度の微分を考える。
\begin{align}
\frac{\partial E(\boldsymbol{w})}{\partial \boldsymbol{w}}&=-\sum_{i=1}^n\frac{\partial E_i}{\partial p_i}・\frac{\partial p_i}{\partial z_i}・\frac{\partial z_i}{\partial \boldsymbol{w}}\\
&=-\sum_{i=1}^n\bigr(\frac{y_i}{p_i}-\frac{1-y_i}{1-p_i} \bigl)p_i(1-p_i)\boldsymbol{x}\\
&=-\sum_{i=1}^n(y_i(1-p_i)-(1-y_i)p_i)\boldsymbol{x}\\
&=-\sum_{i=1}^n(y_i-p_i)\boldsymbol{x}
\end{align}
このとき勾配降下法による重みパラメータの更新は、
\begin{align}
\boldsymbol{w}(k+1)&=\boldsymbol{w}(k)-\eta\frac{\partial E(\boldsymbol{w})}{\partial \boldsymbol{w}}\\
&=\boldsymbol{w}(k)+\eta\sum_{i=1}^n(y_i-p_i)\boldsymbol{x}
\end{align}
と書くことができ、更新されなくなったとき勾配が0になったということとなる。
確率的勾配降下法とは、データを１つずつランダムに選びパラメータを更新
\boldsymbol{w}(k+1)=\boldsymbol{w}(k)+\eta(y_i-p_i)\boldsymbol{x}
性能指標
予測結果と真のラベルクロス集計表を混同行列と呼ぶ。
| 正解 | ラベル | ||
| positive | negetive | ||
| 予測 | positive | 真陽性 | 偽陽性 |
| ラベル | negetive | 偽陰性 | 真陽性 |
-
正解率 正解した数/予測対象全体の数
予測対称クラスに偏りがある場合などに問題がある -
再現率 真陽性/本当に陽性であった数(真陽性+偽陰性)
-
適合率 真陽性/陽性と予測した数(真陽性+偽陽性)
再現率と適合率はトレードオフの関係にるため、これらの調和平均であるF値を使うことがある。
F = 1/(1/再現率 + 1/適合率)
3.2 実装
ダミーデータを用いて動作の確認を行う。
n_sample = 100
harf_n_sample = 50
var = .2
def gen_data(n_sample, harf_n_sample):
x0 = np.random.normal(size=n_sample).reshape(-1, 2) - 1.
x1 = np.random.normal(size=n_sample).reshape(-1, 2) + 1.
x_train = np.concatenate([x0, x1])
y_train = np.concatenate([np.zeros(harf_n_sample), np.ones(harf_n_sample)]).astype(np.int)
return x_train, y_train
def plt_data(x_train, y_train):
plt.scatter(x_train[:, 0], x_train[:, 1], c=y_train, facecolor="none", edgecolor="b", s=50, label="training data")
plt.legend()
#データ作成
x_train, y_train = gen_data(n_sample, harf_n_sample)
#データ表示
plt_data(x_train, y_train)
予測モデルの作成を行い、分類結果を図示する。
データ点以外の領域についても分類結果を出力することでモデルの様子を確認できる。
from sklearn.linear_model import LogisticRegression
xx0, xx1 = np.meshgrid(np.linspace(-5, 5, 100), np.linspace(-5, 5, 100))
xx = np.array([xx0, xx1]).reshape(2, -1).T
model=LogisticRegression(fit_intercept=True)
model.fit(x_train, y_train)
proba = model.predict_proba(xx)
y_pred = (proba > 0.5).astype(np.int)
plt.scatter(x_train[:, 0], x_train[:, 1], c=y_train)
plt.contourf(xx0, xx1, proba[:, 0].reshape(100, 100), alpha=0.2, levels=np.linspace(0, 1, 3))
タイタニックの生存データを用いた予測モデルの構築を行う。
課題としては「30歳の男性の生き残る確率」を求めることである。
import pandas as pd
data_train = pd.read_csv('C:/Users/tanak/study/data/titanic_train.csv')
data_train.head()
目的変数は「Suevived」であり、1が生存、0が死亡である。
今回は、簡単のため説明変数としてAge(年齢)とSex(性別)のみ扱うものとする。
y_train = data_train['Survived'].values
x_train = data_train[['Age', 'Sex']]
欠損値の確認を行う。
x_train.isna().sum()
Ageの欠損値は平均値を使って補完を行うことにする。
x_train.loc[x_train['Age'].isna(), 'Age'] = int(x_train['Age'].mean())
x_train.isna().sum()
Sexが文字列のままであるとモデルの作成が行えないので、maleなら1、femaleなら0とする。
x_train.loc[x_train['Sex']=='male', 'Sex'] = 1
x_train.loc[x_train['Sex']=='female', 'Sex'] = 0
from sklearn.metrics import f1_score, accuracy_score, confusion_matrix
LR = LogisticRegression()
LR.fit(x_train, y_train)
pred_train = LR.predict(x_train)
print('正解率: ', accuracy_score(y_train, pred_train))
print('F1値: ', f1_score(y_train, pred_train))
print('混同行列')
confusion_matrix(y_train, pred_train)
30代男性の生き残る確率を出力する。
LR.predict_proba([[30, 1]])[0][0]
4. 主成分分析
4.1 要約
多変量のデータをより少数のデータに圧縮
分散が最大となる射影軸を探索
学習データを、
\boldsymbol{x_i}=(x_{i1},x_{i2},\cdots,x_{im})
としたとき、平均ベクトルは
\hat{\boldsymbol{x}}=\frac{1}{n}\sum_{i=1}^{n}\boldsymbol{x_i}
と表すことができる。
平均値で中心化を行ったデータ行列を,
\bar{\boldsymbol{X}}=(\boldsymbol{x}_1-\hat{\boldsymbol{x}},\cdots,\boldsymbol{x}_n-\hat{\boldsymbol{x}})
としたときに、分散共分散行列は
\Sigma = Var(\bar{\boldsymbol{X}})=\frac{1}{n}\bar{\boldsymbol{X}}^T\bar{\boldsymbol{X}}
線形変換後のベクトルは、
\boldsymbol{s_j}=(s_{1j},\cdots,s_{nj})=\bar{\boldsymbol{X}}\boldsymbol{a}_j
となる。
線形変換後の分散を計算すると、
Var(\boldsymbol{s_j})=\frac{1}{n}\boldsymbol{s_j^T}\boldsymbol{s_j}=\frac{1}{n}\boldsymbol{a_j^T}\bar{\boldsymbol{X}}^T\bar{\boldsymbol{X}}\boldsymbol{a_j}=\boldsymbol{a_j^T}Var(\bar{\boldsymbol{X}})\boldsymbol{a_j}
以下の制約付き最適化問題を解くことで、線形変換後の分散が最大となる射影軸を求める。
\arg\max_{\boldsymbol{a}\in R^m}\ \boldsymbol{a}_j^TVar(\bar{\boldsymbol{X}})\boldsymbol{a_j}\\
subject\ to\ \boldsymbol{a_j^T}\boldsymbol{a_j}=1
詳細は割愛するが、この解は、
Var(\bar{\boldsymbol{X}})\boldsymbol{a_j}=\lambda_j\boldsymbol{a_j}
となり、射影先の分散は固有値と一致することが確認できる。
Var(\boldsymbol{s_j})=\boldsymbol{a_j^T}Var(\bar{\boldsymbol{X}})\boldsymbol{a_j}=\lambda_j\boldsymbol{a_j^T}\boldsymbol{a_j}=\lambda_j
寄与率
第1～元次元の主成分の分散は、元のデータの分散と一致。
V_{total}=\sum_{i=1}^m\lambda_i
寄与率($c_k$)は、第k成分の分散の全分散に対する割合のことである。
累積寄与率($r_k$)は、第1～第k主成分まで圧縮した際の情報損失量の割合である。
c_k=\frac{\lambda_k}{\sum_{i=1}^m\lambda_i}\\
r_k=\frac{\sum_{i=j}^k\lambda_j}{\sum_{i=1}^m\lambda_i}
4.2 実装
乳がん検査データを用いたロジスティック回帰を行う。
主成分成分分析により、2次元空間上に圧縮し分類できるか確認する。
from sklearn.datasets import load_breast_cancer
from sklearn.preprocessing import StandardScaler
from sklearn.decomposition import PCA
from sklearn.linear_model import LogisticRegression
cancer = load_breast_cancer()
scaler = StandardScaler()
scaler.fit(cancer.data)
y = cancer.target
X_std = scaler.transform(cancer.data)
pca = PCA(n_components=2, whiten=True) # whiten 白色化
lr = LogisticRegression(solver='lbfgs')
X_pca = pca.fit_transform(X_std)
lr.fit(X_pca, y)
xx0, xx1 = np.meshgrid(np.linspace(-5, 6, 100), np.linspace(-5, 6, 100))
xx = np.array([xx0, xx1]).reshape(2, -1).T
proba = lr.predict_proba(xx)
y_pred = (proba > 0.5).astype(np.int)
plt.scatter(X_pca[:, 0], X_pca[:, 1], c=y)
plt.contourf(xx0, xx1, proba[:, 0].reshape(100, 100), alpha=0.2, levels=np.linspace(0, 1, 3))
寄与率を確認すると、2つの主成分で全体の分散のおよそ63%の分散を捉えていることがわかる。
pca.explained_variance_ratio_
5. サポートベクターマシン
5.1 要約
SVMの目的は、マージンの最大化である。
マージンとは、決定境界と、この境界に最も近い訓練データの間の距離として定義される。
この決定境界に最も近い訓練データはサポートベクトルと呼ばれる。
ロジスティック回帰に比べ外れ値の影響を受けにくい。
ここで、決定境界に沿った正と負の超平面を見ていく。
これらの超平面は、次のように表せる。
w_0+\boldsymbol{w}^T\boldsymbol{x}_{pos}=1\\
w_0+\boldsymbol{w}^T\boldsymbol{x}_{neg}=-1
この2つの方程式の差をとる。
\boldsymbol{w}^T(\boldsymbol{x}_{pos}-\boldsymbol{x}_{neg})=2
この式の両辺をベクトルの長さ$|\boldsymbol{w}|$で割る。
\frac{\boldsymbol{w}^T(\boldsymbol{x}_{pos}-\boldsymbol{x}_{neg})}{|\boldsymbol{w}|}=\frac{2}{|\boldsymbol{w}|}
この式の左辺は、正の超平面と負の超平面の距離(すなわちマージン)として解釈できる。
SVMの目的関数を最大化する問題は、$\frac{2}{|\boldsymbol{w}|}$を最大化する問題に帰着する。
実際には、逆数をとって2乗した$\frac{1}{2}|\boldsymbol{w}|^2$を最小化するほうが簡単である。
ただし、これにはデータ点が正しく分類されているという制約がある。
w_0+\boldsymbol{w}^T\boldsymbol{x}^{(i)}\geq1\ (y^{(i)}=1)\\
w_0+\boldsymbol{w}^T\boldsymbol{x}^{(i)}\leq-1\ (y^{(i)}=-1)\\
i=1\cdots N
簡潔に書くと、
y^{(i)}(w_0+\boldsymbol{w}^T\boldsymbol{x}^{(i)})\geq1
となる。
ソフトマージン分類
スラック変数$\xi$を導入することで、適切なペナルティを科した上で、
誤分類が存在する状態のまま最適化問題を収束させることが可能となった。
w_0+\boldsymbol{w}^T\boldsymbol{x}^{(i)}\geq1-\xi^{(i)}\ (y^{(i)}=1)\\
w_0+\boldsymbol{w}^T\boldsymbol{x}^{(i)}\leq-1+\xi^{(i)}\ (y^{(i)}=-1)\\
i=1\cdots N
最小化すべき対象は次のようになる。
\frac{1}{2}|\boldsymbol{w}|^2+C\biggl(\sum_{i}\xi^{(i)} \biggl)
Cの値が大きい場合は、誤分類のペナルティが大きく、Cの値が小さいときはペナルティが小さい。
これにより、Cを使ってマージンの幅を制御できる。
カーネルSVM
線形分離不能なデータの処理をする方法に、カーネル法がある。
射影関数$\phi$を使ってそれらの組み合わせを高次元へ射影し、線形分離できるようにする。
例えば、
\phi(x_1,x_2)=(z_1,z_2,z_3)=(x_1,x_2,x_1^2+x_2^2)
のような変換である。
ただし、新しい特徴量を生成する計算コストが非常に高いという問題がある。
ここで、カーネルトリックと呼ばれる手法があり、カーネル関数を定義して、データ間の距離を直接計算する。
K(\boldsymbol{x}^{(i)}, \boldsymbol{x}^{(j)})=\phi(\boldsymbol{x}^{(i)})^T\phi(\boldsymbol{x}^{(j)})
最もよく使われるカーネルの1つに**動径基底関数カーネル(ガウスカーネル)**がある。
K(\boldsymbol{x}^{(i)}, \boldsymbol{x}^{(j)})=\exp\biggl(-\frac{|\boldsymbol{x}^{(i)}-\boldsymbol{x}^{(j)}|^2}{2\sigma^2} \biggr)
これは、次のように簡略化される。
K(\boldsymbol{x}^{(i)}, \boldsymbol{x}^{(j)})=\exp\bigl(-\gamma|\boldsymbol{x}^{(i)}-\boldsymbol{x}^{(j)}|^2 \bigr)
ここで、$\gamma=\frac{1}{2\sigma^2}$はハイパーパラメータである。
5.2 実装
numpyによる実装とschikit-learnによる実装の両方を乗せる。
線形分離可能
def gen_data():
x0 = np.random.normal(size=50).reshape(-1, 2) - 2.
x1 = np.random.normal(size=50).reshape(-1, 2) + 2.
X_train = np.concatenate([x0, x1])
ys_train = np.concatenate([np.zeros(25), np.ones(25)]).astype(np.int)
return X_train, ys_train
X_train, ys_train = gen_data()
plt.scatter(X_train[:, 0], X_train[:, 1], c=ys_train)
t = np.where(ys_train == 1.0, 1.0, -1.0)
n_samples = len(X_train)
# 線形カーネル
K = X_train.dot(X_train.T)
eta1 = 0.01
eta2 = 0.001
n_iter = 500
H = np.outer(t, t) * K
a = np.ones(n_samples)
for _ in range(n_iter):
grad = 1 - H.dot(a)
a += eta1 * grad
a -= eta2 * a.dot(t) * t
a = np.where(a > 0, a, 0)
index = a > 1e-6
support_vectors = X_train[index]
support_vector_t = t[index]
support_vector_a = a[index]
term2 = K[index][:, index].dot(support_vector_a * support_vector_t)
b = (support_vector_t - term2).mean()
xx0, xx1 = np.meshgrid(np.linspace(-5, 5, 100), np.linspace(-5, 5, 100))
xx = np.array([xx0, xx1]).reshape(2, -1).T
X_test = xx
y_project = np.ones(len(X_test)) * b
for i in range(len(X_test)):
for a, sv_t, sv in zip(support_vector_a, support_vector_t, support_vectors):
y_project[i] += a * sv_t * sv.dot(X_test[i])
y_pred = np.sign(y_project)
# 訓練データを可視化
plt.scatter(X_train[:, 0], X_train[:, 1], c=ys_train)
# サポートベクトルを可視化
plt.scatter(support_vectors[:, 0], support_vectors[:, 1],
s=100, facecolors='none', edgecolors='k')
# 領域を可視化
#plt.contourf(xx0, xx1, y_pred.reshape(100, 100), alpha=0.2, levels=np.linspace(0, 1, 3))
# マージンと決定境界を可視化
plt.contour(xx0, xx1, y_project.reshape(100, 100), colors='k',
levels=[-1, 0, 1], alpha=0.5, linestyles=['--', '-', '--'])
# マージンと決定境界を可視化
plt.quiver(0, 0, 0.1, 0.35, width=0.01, scale=1, color='pink')
from sklearn.svm import SVC
svc = SVC(kernel='linear', degree=2, probability=True) # Cはペナルティ
svc.fit(X_train, ys_train)
xx0, xx1 = np.meshgrid(np.linspace(-5, 5, 100), np.linspace(-5, 5, 100))
xx = np.array([xx0, xx1]).reshape(2, -1).T
X_test = xx
y_pred = svc.predict(X_test)
support_vectors = svc.support_vectors_
# 訓練データを可視化
plt.scatter(X_train[:, 0], X_train[:, 1], c=ys_train)
# サポートベクトルを可視化
plt.scatter(support_vectors[:, 0], support_vectors[:, 1],
s=100, facecolors='none', edgecolors='k')
# 領域を可視化
#plt.contourf(xx0, xx1, y_pred.reshape(100, 100), alpha=0.2, levels=np.linspace(0, 1, 3))
# マージンと決定境界を可視化
plt.contour(xx0, xx1, svc.decision_function(X_test).reshape(100, 100), colors='k',
levels=[-1, 0, 1], alpha=0.5, linestyles=['--', '-', '--'])
# マージンと決定境界を可視化
plt.quiver(0, 0, 0.1, 0.35, width=0.01, scale=1, color='pink')
ソフトマージンSVM
x0 = np.random.normal(size=50).reshape(-1, 2) - 1.
x1 = np.random.normal(size=50).reshape(-1, 2) + 1.
x_train = np.concatenate([x0, x1])
y_train = np.concatenate([np.zeros(25), np.ones(25)]).astype(np.int)
plt.scatter(x_train[:, 0], x_train[:, 1], c=y_train)
X_train = x_train
t = np.where(y_train == 1.0, 1.0, -1.0)
n_samples = len(X_train)
# 線形カーネル
K = X_train.dot(X_train.T)
C = 1
eta1 = 0.01
eta2 = 0.001
n_iter = 1000
H = np.outer(t, t) * K
a = np.ones(n_samples)
for _ in range(n_iter):
grad = 1 - H.dot(a)
a += eta1 * grad
a -= eta2 * a.dot(t) * t
# ここが変わる
a = np.clip(a, 0, C)
index = a > 1e-8
support_vectors = X_train[index]
support_vector_t = t[index]
support_vector_a = a[index]
term2 = K[index][:, index].dot(support_vector_a * support_vector_t)
b = (support_vector_t - term2).mean()
xx0, xx1 = np.meshgrid(np.linspace(-5, 5, 100), np.linspace(-5, 5, 100))
xx = np.array([xx0, xx1]).reshape(2, -1).T
X_test = xx
y_project = np.ones(len(X_test)) * b
for i in range(len(X_test)):
for a, sv_t, sv in zip(support_vector_a, support_vector_t, support_vectors):
y_project[i] += a * sv_t * sv.dot(X_test[i])
y_pred = np.sign(y_project)
# 訓練データを可視化
plt.scatter(x_train[:, 0], x_train[:, 1], c=y_train)
# サポートベクトルを可視化
plt.scatter(support_vectors[:, 0], support_vectors[:, 1],
s=100, facecolors='none', edgecolors='k')
# 領域を可視化
plt.contourf(xx0, xx1, y_pred.reshape(100, 100), alpha=0.2, levels=np.linspace(0, 1, 3))
# マージンと決定境界を可視化
plt.contour(xx0, xx1, y_project.reshape(100, 100), colors='k',
levels=[-1, 0, 1], alpha=0.5, linestyles=['--', '-', '--'])
svc = SVC(C=1, kernel='linear', degree=2, probability=True) # Cはペナルティ
svc.fit(X_train, ys_train)
xx0, xx1 = np.meshgrid(np.linspace(-5, 5, 100), np.linspace(-5, 5, 100))
xx = np.array([xx0, xx1]).reshape(2, -1).T
X_test = xx
y_pred = svc.predict(X_test)
support_vectors = svc.support_vectors_
# 訓練データを可視化
plt.scatter(X_train[:, 0], X_train[:, 1], c=ys_train)
# サポートベクトルを可視化
plt.scatter(support_vectors[:, 0], support_vectors[:, 1],
s=100, facecolors='none', edgecolors='k')
# 領域を可視化
#plt.contourf(xx0, xx1, y_pred.reshape(100, 100), alpha=0.2, levels=np.linspace(0, 1, 3))
# マージンと決定境界を可視化
plt.contour(xx0, xx1, svc.decision_function(X_test).reshape(100, 100), colors='k',
levels=[-1, 0, 1], alpha=0.5, linestyles=['--', '-', '--'])
カーネルSVM
factor = .2
n_samples = 50
linspace = np.linspace(0, 2 * np.pi, n_samples // 2 + 1)[:-1]
outer_circ_x = np.cos(linspace)
outer_circ_y = np.sin(linspace)
inner_circ_x = outer_circ_x * factor
inner_circ_y = outer_circ_y * factor
X = np.vstack((np.append(outer_circ_x, inner_circ_x),
np.append(outer_circ_y, inner_circ_y))).T
y = np.hstack([np.zeros(n_samples // 2, dtype=np.intp),
np.ones(n_samples // 2, dtype=np.intp)])
X += np.random.normal(scale=0.15, size=X.shape)
x_train = X
y_train = y
plt.scatter(x_train[:,0], x_train[:,1], c=y_train)
def rbf(u, v):
sigma = 0.8
return np.exp(-0.5 * ((u - v)**2).sum() / sigma**2)
X_train = x_train
t = np.where(y_train == 1.0, 1.0, -1.0)
n_samples = len(X_train)
# RBFカーネル
K = np.zeros((n_samples, n_samples))
for i in range(n_samples):
for j in range(n_samples):
K[i, j] = rbf(X_train[i], X_train[j])
eta1 = 0.01
eta2 = 0.001
n_iter = 5000
H = np.outer(t, t) * K
a = np.ones(n_samples)
for _ in range(n_iter):
grad = 1 - H.dot(a)
a += eta1 * grad
a -= eta2 * a.dot(t) * t
a = np.where(a > 0, a, 0)
index = a > 1e-6
support_vectors = X_train[index]
support_vector_t = t[index]
support_vector_a = a[index]
term2 = K[index][:, index].dot(support_vector_a * support_vector_t)
b = (support_vector_t - term2).mean()
xx0, xx1 = np.meshgrid(np.linspace(-1.5, 1.5, 100), np.linspace(-1.5, 1.5, 100))
xx = np.array([xx0, xx1]).reshape(2, -1).T
X_test = xx
y_project = np.ones(len(X_test)) * b
for i in range(len(X_test)):
for a, sv_t, sv in zip(support_vector_a, support_vector_t, support_vectors):
y_project[i] += a * sv_t * rbf(X_test[i], sv)
y_pred = np.sign(y_project)
# 訓練データを可視化
plt.scatter(x_train[:, 0], x_train[:, 1], c=y_train)
# サポートベクトルを可視化
plt.scatter(support_vectors[:, 0], support_vectors[:, 1],
s=100, facecolors='none', edgecolors='k')
# 領域を可視化
plt.contourf(xx0, xx1, y_pred.reshape(100, 100), alpha=0.2, levels=np.linspace(0, 1, 3))
# マージンと決定境界を可視化
plt.contour(xx0, xx1, y_project.reshape(100, 100), colors='k',
levels=[-1, 0, 1], alpha=0.5, linestyles=['--', '-', '--'])
svc = SVC(C=1, kernel='rbf', degree=2, probability=True) # Cはペナルティ
svc.fit(X_train, ys_train)
xx0, xx1 = np.meshgrid(np.linspace(-1.5, 1.5, 100), np.linspace(-1.5, 1.5, 100))
xx = np.array([xx0, xx1]).reshape(2, -1).T
X_test = xx
y_pred = svc.predict(X_test)
support_vectors = svc.support_vectors_
# 訓練データを可視化
plt.scatter(X_train[:, 0], X_train[:, 1], c=ys_train)
# サポートベクトルを可視化
plt.scatter(support_vectors[:, 0], support_vectors[:, 1],
s=100, facecolors='none', edgecolors='k')
# 領域を可視化
#plt.contourf(xx0, xx1, y_pred.reshape(100, 100), alpha=0.2, levels=np.linspace(0, 1, 3))
# マージンと決定境界を可視化
plt.contour(xx0, xx1, svc.decision_function(X_test).reshape(100, 100), colors='k',
levels=[-1, 0, 1], alpha=0.5, linestyles=['--', '-', '--'])
