0
1

Delete article

Deleted articles cannot be recovered.

Draft of this article would be also deleted.

Are you sure you want to delete this article?

More than 3 years have passed since last update.

ラビット・チャレンジ(機械学習レポート)

Last updated at Posted at 2021-07-17

1.線形回帰モデル

回帰問題とは、ある入力(離散あるいは連続値)から出力(連続値)を予測する問題

線形回帰 非線形回帰
直線で予測 曲線で予測

線形回帰モデル

・回帰問題を解くための機械学習モデルのひとつ
・教師あり学習
・入力とm次元パラメータの線型結合を出力するモデル

※慣例として予測値にはハット($\hat{ }$)をつける(正解データと区別するため)

説明変数:入力(m次元のベクトル(m=1の場合はスカラ))

\boldsymbol{x}=(x_1, x_2, \cdots, x_m)^T\in\mathbb{R}^m

目的変数:出力(スカラー値)

y\in\mathbb{R}^1
\begin{aligned}\hat{y}&=\boldsymbol{w}^T\boldsymbol{x}+w_0\\
&=
\begin{pmatrix}
w_1 & w_2 & \cdots & w_m
\end{pmatrix}
\begin{pmatrix}
x_1\\
x_2\\
\cdots\\
x_m
\end{pmatrix}+w_0\\
&=\{w_1x_1+w_2x_2+\cdots+w_mx_m\}+w_0\\
&=\sum_{j=1}^{m}w_jx_j+w_0
\end{aligned}

・入力ベクトルと未知のパラメータの各要素を掛け算し、足し合わせたもの
・入力ベクトルとの線型結合に加え、切片も足し合わせる
・入力ベクトルは多次元でも、出力はスカラ(1次元)となる

回帰係数

\boldsymbol{w}=(w_1, w_2, \cdots, w_m)^T\in\mathbb{R}^m

・モデルに含まれる推定すべき未知のパラメータ
・特徴量が予測値に対してどのように影響を与えるかを決定する重みの集合
 (重みが大きければ、その特徴量は予測に大きな影響力を持つ)
   - 重みが正:その特徴量を増加→予測値が増加
   - 重みが0:その特徴量は予測に全く影響しない
   - 重みが負:その特徴量を増加→予測値が減少
・切片:y軸との交点を表す

教師データ

~\{(\boldsymbol{x}_i, y_i);i=1,\cdots, n\}~

線形単回帰モデル

・説明変数が1次元(m=1)の場合、単回帰モデルという
・データは回帰直線に誤差が加わり観測されていると仮定

モデル数式と幾何学的意味

\large y=w_0+w_1x_1+\epsilon\\
\small y\tiny~(目的変数),~~~
\small w_0\tiny~(切片),~~~
\small w_1\tiny~(回帰係数),~~~
\small x_1\tiny~(説明変数(特徴量)),~~~
\small \epsilon\tiny~(誤差)\qquad\qquad\qquad

<参考動画>

線形重回帰モデル

・説明変数が多次元の場合(m>1)の場合、重回帰モデルという
・データは回帰平面(回帰超平面)に誤差が加わり観測されていると仮定

モデル数式と幾何学的意味

\large y=w_0+w_1x_1+w_2x_2+\cdots+w_{n-1}x_{n-1}+\epsilon\\
\small y\tiny~(目的変数),~~~
\small w_0\tiny~(切片),~~~
\small w_1,w_2\tiny~(回帰係数),~~~
\small x_1,x_2\tiny~(説明変数(特徴量)),~~~
\small \epsilon\tiny~(誤差)\qquad\qquad
\begin{align}
y&=w_0+w_1x_1+w_2x_2+\cdots+w_{n-1}x_{n-1}+\epsilon\\
&=w_0+\sum_{i=1}^{n-1}w_ix_i\quad(w_0=w_0\times1=w_0x_0)\\
&=\sum_{i=0}^{n-1}w_ix_i\quad ,{\rm where}~~x_0\equiv1\\
&=\boldsymbol{w}^T\boldsymbol{x}\quad ,{\rm where}~~\boldsymbol{w}=
\begin{pmatrix}
w_0\\
w_1\\
\vdots\\
w_{n-1}
\end{pmatrix}\\
&~~~\qquad\qquad\qquad\boldsymbol{w}^T=
\begin{pmatrix}
w_0 & w_1 \cdots & w_{n-1}
\end{pmatrix}\\
&~~~\qquad\qquad\qquad\boldsymbol{x}=
\begin{pmatrix}
x_0\\
x_1\\
\vdots\\
x_{n-1}
\end{pmatrix}
\end{align}

連立方程式と行列表現

連立方程式

<線形単回帰モデル>

\begin{equation}
\left\{ \,
    \begin{aligned}
    & y_1=w_0+w_1x_1+\epsilon_1 \\
    & y_2=w_0+w_1x_2+\epsilon_2 \\
    & \quad\vdots \\
    & y_n=w_0+w_1x_n+\epsilon_n
    \end{aligned}
\right.
\end{equation}

<線形重回帰モデル>

\begin{equation}
\left\{ \,
    \begin{aligned}
    & y_1=w_0+w_1x_{11}+w_2x_{12}+\cdots+w_mx_{1m}+\epsilon_1 \\
    & y_1=w_0+w_1x_{21}+w_2x_{22}+\cdots+w_mx_{2m}+\epsilon_2 \\
    & \quad\vdots \\
    & y_1=w_0+w_1x_{n1}+w_2x_{n2}+\cdots+w_mx_{nm}+\epsilon_n
    \end{aligned}
\right.
\end{equation}

※この式の$m=1$のときが線形単回帰モデルである

行列表現

線形単回帰モデルも線形重回帰モデルも、行列表現をすると以下となる。
この式の行列$X$を__計画行列__という

\Large \boldsymbol{y}=X\boldsymbol{w}+\boldsymbol{\epsilon} \\

<線形単回帰モデル>
行列の要素で具体的に表現すると以下のようになり、線形単回帰モデルの連立方程式と一致することがわかる

\underset{n\times 1}{\begin{pmatrix}
 y_1 \\
 y_2 \\
 \vdots \\
 y_n
\end{pmatrix}}
=\underset{n\times 2}{\begin{pmatrix}
 1 & x_1 \\
 1 & x_2 \\
 \vdots & \vdots \\
 1 & x_n
\end{pmatrix}} 
\underset{2\times 1}{\begin{pmatrix}
 ~ \\
 w_0 \\
 w_1 \\
 ~
\end{pmatrix}}
+\underset{n\times 1}{\begin{pmatrix}
 \epsilon_1 \\
 \epsilon_2 \\
 \vdots \\
 \epsilon_n
\end{pmatrix}}

<線形重回帰モデル>
行列の要素で具体的に表現すると以下のようになり、線形重回帰モデルの連立方程式と一致することがわかる

\underset{n\times 1}{\begin{pmatrix}
 y_1 \\
 y_2 \\
 \vdots \\
 y_n
\end{pmatrix}}
=\underset{n\times (m+1)}{\begin{pmatrix}
 1 & x_{11} & x_{12} & \cdots & x_{1m} \\
 1 & x_{21} & x_{22} & \cdots & x_{2m} \\
 \vdots & \vdots & \vdots & & \vdots \\
 1 & x_{n1} & x_{n2} & \cdots & x_{nm} 
\end{pmatrix}} 
\underset{(m+1)\times 1}{\begin{pmatrix}
 w_0 \\
 w_1 \\
 \vdots \\
 w_m
\end{pmatrix}}
+\underset{n\times 1}{\begin{pmatrix}
 \epsilon_1 \\
 \epsilon_2 \\
 \vdots \\
 \epsilon_n
\end{pmatrix}}

データの分割とモデルの汎化性能測定

データの分割

データの分割はモデルの汎化性能(Generalization)を測定するために行うものであり、単にデータへのあてはまりの良さというわけではなく、未知のデータに対してどれくらい精度が高いかを測るために行うものである

学習用データ:機械学習モデルの学習に利用するデータ
検証用データ:学習済みモデルの精度を検証するためのデータ

モデルの汎化性能測定

  1. すべてのデータ $\{ (\boldsymbol x_i, y_i);i=1, \cdots, n\}~$ を学習用データと検証用データに分割
  2. 学習用データ $\{ (\boldsymbol x_i^{(train)} , y_i^{(train)});i=1, \cdots, n_{train}\}$ を利用して推定パラメータを学習
  3. 検証用データ $\{ (\boldsymbol x_i^{(test)}, y_i^{(test)});i=1, \cdots, n_{test}\}$ で検証

学習誤差

MSE_{train}=\frac{1}{n_{train}} \sum_{i=1}^{n_{train}}\bigl(\hat{y}_i^{(train)}-y_i^{(train)}\bigr)^2~

検証誤差

MSE_{test}=\frac{1}{n_{test}} \sum_{i=1}^{n_{test}}\bigl(\hat{y}_i^{(test)}-y_i^{(test)}\bigr)^2

線形回帰モデルのパラメータ推定

平均二乗誤差(残差平方和)

データ $y_i^{(train)}$ とモデルの出力 $\hat{y}_i^{(train)}$ の二乗誤差の和

  • 小さいほど回帰直線(回帰超平面)とデータの距離が近い
  • パラメータのみに依存する関数(データ $\boldsymbol x_i^{(train)} , y_i^{(train)}$ は既知、パラメータ $\boldsymbol{w}$ のみ未知)
    • $\boldsymbol{w}$だけを変数とみた$\boldsymbol{w}$の関数ということ
\begin{align}
MSE_{train}&=\frac{1}{n_{train}} \sum_{i=1}^{n_{train}}\bigl(\hat{y}_i^{(train)}-y_i^{(train)}\bigr)^2~\\
&=\frac{1}{n_{train}} \sum_{i=1}^{n_{train}}\bigl(\boldsymbol{w}^T\boldsymbol{x}_i^{(train)}-y_i^{(train)}\bigr)^2
\end{align}
  • $MSE$ は、$J(\boldsymbol{w})$ とか $E(\boldsymbol{w})$ と書かれることも多い
  • 二乗損失は一般に外れ値に弱い

最小二乗法

学習データの平均二乗誤差を最小とするパラメータ(回帰係数)を探索
(MSEを最小にするような $\hat{\boldsymbol{w}}$ (m+1次元)を求める)

\large \hat{\boldsymbol{w}}=\arg \underset{\boldsymbol{w} \in \mathbb{R}^{m+1}}{\min} MSE_{train}

学習データの平均二乗誤差の最小化は、その勾配が0になる点を求めればよい
(MSEを $\boldsymbol{w}$ で偏微分したものが0となる $\boldsymbol{w}$ の点を求める)

\large\frac{\partial}{\partial \boldsymbol{w}}MSE_{train}=0
__回帰係数__ $\hat{\boldsymbol{w}}$(ここをクリックして計算式を表示)
\begin{align}
\Longrightarrow\quad &\frac{\partial}{\partial \boldsymbol{w}} \biggl\{\frac{1}{n_{train}}  \sum_{i=1}^{n_{train}}\bigl(\hat{y}_i^{(train)}-y_i^{(train)}\bigr)^2 \biggr\} =0 \\
\Longrightarrow\quad &\frac{\partial}{\partial \boldsymbol{w}} \biggl\{\frac{1}{n_{train}} \sum_{i=1}^{n_{train}}\bigl(\boldsymbol{w}^T\boldsymbol{x}_i^{(train)}-y_i^{(train)}\bigr)^2\biggr\} =0 \\
\Longrightarrow\quad &\frac{\partial}{\partial \boldsymbol{w}} \biggl\{\frac{1}{n_{train}} \sum_{i=1}^{n_{train}}\bigl({\boldsymbol{x}_i^{(train)}}^T\boldsymbol{w}-y_i^{(train)}\bigr)^2\biggr\} =0 \\
\Longrightarrow\quad &\frac{\partial}{\partial \boldsymbol{w}} \biggl\{\frac{1}{n_{train}} (X\boldsymbol{w}-\boldsymbol{y})^T(X\boldsymbol{w}-\boldsymbol{y})\biggr\} =0 \\
\Longrightarrow\quad &\frac{1}{n_{train}} \frac{\partial}{\partial \boldsymbol{w}} \bigl\{ (X\boldsymbol{w}-\boldsymbol{y})^T(X\boldsymbol{w}-\boldsymbol{y})\bigr\} =0 \\
\Longrightarrow\quad &\frac{1}{n_{train}} \frac{\partial}{\partial \boldsymbol{w}} \bigl\{ (\boldsymbol{w}^TX^T-\boldsymbol{y}^T)(X\boldsymbol{w}-\boldsymbol{y})\bigr\} =0 \\
\Longrightarrow\quad &\frac{1}{n_{train}} \frac{\partial}{\partial \boldsymbol{w}} \bigl( \boldsymbol{w}^TX^TX\boldsymbol{w}-\boldsymbol{w}^TX^T\boldsymbol{y}-\boldsymbol{y}^TX\boldsymbol{w}+\boldsymbol{y}^T\boldsymbol{y}\bigr) =0 \\
\Longrightarrow\quad &\frac{1}{n_{train}} \frac{\partial}{\partial \boldsymbol{w}} \bigl( \boldsymbol{w}^TX^TX\boldsymbol{w}-2\boldsymbol{w}^TX^T\boldsymbol{y}+\boldsymbol{y}^T\boldsymbol{y}\bigr) =0 \\
\Longrightarrow\quad &\frac{1}{n_{train}} \bigl( 2X^TX\boldsymbol{w}-2X^T\boldsymbol{y}\bigr) =0 \\

&\quad \because \frac{\partial}{\partial \boldsymbol{w}} \boldsymbol{w}^T A \boldsymbol{w} = (A+A^T)\boldsymbol{w}=2A\boldsymbol{w}~(A:対称行列) \\

&~\quad\quad \frac{\partial}{\partial \boldsymbol{w}} \boldsymbol{w}^T\boldsymbol{x}=\boldsymbol{w} \\

\Longrightarrow\quad &2X^TX\boldsymbol{w}-2X^T\boldsymbol{y}=0 \\
\Longrightarrow\quad &2X^TX\boldsymbol{w}=2X^T\boldsymbol{y} \\
\Longrightarrow\quad &X^TX\boldsymbol{w}=X^T\boldsymbol{y} \\
\Longrightarrow\quad &(X^TX)^{-1}(X^TX)\boldsymbol{w}=(X^TX)^{-1}X^T\boldsymbol{y} \\
\Longrightarrow\quad &\boldsymbol{w}=(X^TX)^{-1}X^T\boldsymbol{y} \\
\end{align}~
\hat{\boldsymbol{w}}=\bigl({X^{(train)}}^TX^{(train)}\bigr)^{-1}{X^{(train)}}^T\boldsymbol{y}^{(train)}

これより、予測したい $n_*$ 個の新たな入力点を $X_*$ とすると__予測値__ $\hat{\boldsymbol{y}}$ は以下となる

\begin{align}
\hat{\boldsymbol{y}}&=\underset{n_*~\times~{(m+1)}}{X_*}~\underset{(m+1)\times1}{\hat{\boldsymbol{w}}} \\
&=X_* \bigl({X^{(train)}}^TX^{(train)}\bigr)^{-1}{X^{(train)}}^T\boldsymbol{y}^{(train)}
\end{align}
X_*=
\begin{pmatrix}
1 & x^*_{11} & \cdots & x^*_{1m} \\
1 & x^*_{21} & \cdots & x^*_{2m} \\
\vdots & \vdots &  & \vdots \\
1 & x^*_{n_*1} & \cdots & x^*_{n_*m} \\
\end{pmatrix}\qquad\qquad\quad

最尤法による回帰係数の推定

  • 誤差を正規分布に従う確率変数を仮定した尤度関数の最大化を利用した推定も可能
  • 線形回帰の場合には、最尤法による解は最小二乗法の解と一致

<参考動画>

実装演習

'skl_regression.ipynb' を使って実装演習を実施

'Boston house prices dataset' は、13 個の説明変数(部屋数 RM , 犯罪率 CRIM など)と 1 個の目的変数(住宅価格 target)で構成されている。

線形単回帰分析

単回帰分析なので説明変数として部屋数RMのみを使って単回帰分析を実施し、目的変数の住宅価格(PRICEとする)を推定する。

部屋数を1〜7部屋で住宅価格の予測値を出してみた。
スクリーンショット 2021-07-13 14.20.33.png
部屋数が1〜3部屋の住宅価格がマイナスの値になっており、明らかにおかしい。
これは教師データとなるデータセットに部屋数が少ない住宅のデータが全データに比べて著しく少ないもしくは存在しないのではないかと推測される。

重回帰分析(2変数)

説明変数として部屋数RMと犯罪率CRIMの2つを使って重回帰分析を実施し、目的変数の住宅価格(PRICEとする)を推定する。

部屋数を6部屋に固定し、犯罪率を 0.0 〜 1.0 で住宅価格の予測値を出してみた。
スクリーンショット 2021-07-13 15.34.20.png

犯罪率を変えても住宅価格にはあまり影響しないことがわかる。
次に単回帰分析で確認したように、部屋数を1〜6部屋で住宅価格の予測値を出してみた。

スクリーンショット 2021-07-13 15.34.31.png

単回帰分析と同様に、部屋数が1〜3部屋の住宅価格がマイナスの値になってしまった。

回帰係数と切片の値を確認

単回帰及び重回帰の怪奇係数と切片を確認する。

スクリーンショット 2021-07-13 15.11.42.png

この回帰係数より、部屋数が増えると住宅価格が上がり、犯罪率が高くなると住宅価格が下がるという傾向がわかる。
しかし、単回帰/重回帰いずれも住宅価格というマイナスになるはずのない予測値(回帰直線、回帰超平面)の切片がマイナスの値を取っていることから、直線的な線形回帰には適さないデータである可能性があるということが考えられる。または、モデルを作るときに使用した教師データの説明変数の範囲外の予測をしようとしている可能性がある(外挿)。

2.非線形回帰モデル

基底展開法

  • 回帰関数として、基底関数と呼ばれる既知の非線形関数 $\phi(\boldsymbol{x})$ とパラメータベクトル $\boldsymbol{w}$ の線形結合を使用
y_i=w_0+\sum_{i=1}^m w_j\phi_j(\boldsymbol{x}_i)+\epsilon_i
  • 未知パラメータの推定:線形回帰モデルと同様に最小二乗法や最尤法を使用

よく使われる基底関数と非線形回帰

1次元の基底関数

|多項式関数|ガウス型基底関数|
|:-:|:-:|:-:|
|$\quad\quad~~~\phi_j(x)=x^j\quad\quad ~~~ $|$\phi_i(x)=\exp \biggl\{ -\frac{(x-\mu_i)^2}{2h_j} \biggr\}~$|
|||

2次元の基底関数

2次元ガウス型基底関数
$\phi_j(\boldsymbol{x})=\exp \biggl\{ -\frac{(\boldsymbol{x}-\boldsymbol{\mu}_j)^T(\boldsymbol{x}-\boldsymbol{\mu}_j)}{2h_j} \biggr\}~$

これら以外にも、スプライン関数/Bスプライン関数などがある

非線形回帰モデル

説明変数

\boldsymbol{x}_i = (x_{i1}, x_{i2}, \cdots, x_{im}) \in \mathbb{R}^m\\

非線形関数ベクトル(k次元の特徴ベクトル)

\boldsymbol{\phi}(\boldsymbol{x}_i) = (\phi_1(\boldsymbol{x}_i), \phi_2(\boldsymbol{x}_i), \cdots, \phi_k(\boldsymbol{x}_i))^T \in \mathbb{R}^k\\

非線形関数の計画行列(k次元のパラメータベクトルがn個)

\Phi^{(train)} = (\boldsymbol{\phi}(\boldsymbol{x}_1), \boldsymbol{\phi}(\boldsymbol{x}_2), \cdots, \boldsymbol{\phi}(\boldsymbol{x}_n))^T \in \mathbb{R}^{n\times k}\\

最尤法による予測値

\boldsymbol{\hat{y}}=\Phi(\Phi^{(train)T}\Phi^{(train)})^{-1}\Phi^{(train)T}\boldsymbol{y}^{(train)}\\

未学習と過学習

  • 未学習(underfitting)とは、学習データに対して、十分小さな誤差が得られないモデル
    • 対策 → モデルの表現力が低いため、表現力の高いモデルを利用する
  • 過学習(overfitting)とは、小さな誤差は得られたけど、テスト集合誤差との差が大きいモデル
    • 対策1 → 学習データの数を増やす
    • 対策2 → __不要な基底関数(変数)を削除__して表現力を抑止
    • 対策3 → __正則化法を利用__して表現力を抑止
    • 上記の対策1と2は、モデルの複雑さを調整する方法
      スクリーンショット 2021-07-13 18.51.50.png

過学習対策

① 学習データの数を増やす

サンプル数 100 サンプル数 10000
image.png

正則化パラメータを 0 、基底関数を 50 で固定し、サンプル数のみを変えて比較
サンプル数を増やすことで過学習が回避できたことがわかる

② 不要な基底関数を削除

  • 基底関数の数、位置やバンド幅によりモデルの複雑さが変化
    • 解きたい問題に対して多くの基底関数を用意してしまうと過学習の問題がおこるため、適切
      な基底関数を用意(CVなどで選択)
基底関数 2 基底関数 10 基底関数 50

正則化パラメータを 0 、パラメータ数を 100 で固定し、基底関数のみを変えて比較
この例では、基底関数を 10 としたときに、適切なモデルに近づいていることがわかる

③ 正則化法(罰則化法)

「モデルの複雑さに伴って、その値が大きくなる__正則化項(罰則項)__を課した関数」を最小化

  • 正則化項(罰則項)$R(\boldsymbol{w})$
    形状によっていくつもの種類があり、それぞれ推定量の性質が異なる

  • 正則化(平滑化)パラメータ $\gamma$
    モデルの曲線のなめらかさを調節▶適切に決める必要あり

\large S_\gamma = \underset{MSE}{\underline{(\boldsymbol{y}-\overset{n \times k}{\Phi}\boldsymbol{w})^T(\boldsymbol{y}-\overset{n \times k}{\Phi}\boldsymbol{w})}}+\gamma \underset{正則化項}{\underline{R(\boldsymbol{w})}}\qquad\small(\gamma>0)

基底関数の数(k)が増加すると、パラメータが増加し、残差は減少(モデルが複雑化)

正則化パラメータ 0 正則化パラメータ 0.1

サンプル数を 100 、基底関数を 50 で固定し、正則化パラメータのみを変えて比較
正則化パラメータを 0.1 として正則化項を生かすことにより、過学習が回避できたことがわかる

Ridge推定量(L1正則化)
正則化項 $R(\boldsymbol{w})$ に L2 ノルムを使うのが Ridge 回帰

E(\boldsymbol{w})=(\boldsymbol{y}-{\Phi}\boldsymbol{w})^T(\boldsymbol{y}-\Phi\boldsymbol{w})+\gamma\boldsymbol{w}^T\boldsymbol{w}

Lasso推定量(L2正則化)
正則化項 $R(\boldsymbol{w})$ に L1 ノルムを使うのが Lasso 回帰

E(\boldsymbol{w})=(\boldsymbol{y}-{\Phi}\boldsymbol{w})^T(\boldsymbol{y}-\Phi\boldsymbol{w})+\gamma\sum_{j=1}^d\|\boldsymbol{w}_j\|_1

スクリーンショット 2021-07-13 23.27.17.png

データの分割とモデルの汎化性能測定

汎化性能
学習に使用した入力だけでなく、これまで見たことのない新たな入力に対する予測性能

  • (学習誤差ではなく)汎化誤差(テスト誤差)が小さいものが良い性能を持ったモデル
  • 汎化誤差は通常、学習データとは別に収集された検証データでの性能を測ることで推定

ホールドアウト法

有限のデータを学習用とテスト用の2つに分割し、「予測精度」や「誤り率」を推定する為に使用

  • 学習用を多くすれば学習精度は良くなるが、テスト用が減り性能評価の精度は悪くなる。
  • 逆にテスト用を多くすれば、学習用が減少するので学習そのものの精度が悪くなることになる。
  • __手元にデータが大量にある場合を除いて、良い性能評価を与えないという欠点__がある。

基底展開法に基づく非線形回帰モデルでは、基底関数の数、位置、バンド幅の値とチューニングパラメータをホールドアウト値を小さくするモデルで決定する

クロスバリデーション(交差検証)

データを学習用と評価用に分割(5分割の例)
スクリーンショット 2021-07-14 0.20.10.png

検証データで各モデルの精度を計測
スクリーンショット 2021-07-14 0.20.54.png

CV値の最もよいモデルを採用
スクリーンショット 2021-07-14 0.22.56.png

グリッドサーチ

  • 全てのチューニングパラメータの組み合わせで評価値を算出
  • 最も良い評価値を持つチューニングパラメータを持つ組み合わせを、「いいモデルのパラメータ」として採用

<参考動画>

実装演習

以下の真の関数にノイズを加えてデータを生成

z = 1-48x+218x^2-315x^3+145x^4

スクリーンショット 2021-07-14 3.24.38.png

線形回帰モデル

sklearn.linear_model.LinearRegression
スクリーンショット 2021-07-14 10.50.26.png

決定係数が 0.238 程度であり、線形回帰のグラフなので直線となり、グラフでも真の関数とはかけ離れていることがわかる。

カーネルリッジ回帰(KRR)

sklearn.kernel_ridge.KernelRidge
スクリーンショット 2021-07-14 10.50.38.png

決定係数が 0.829 程度となり、グラフでも真の関数に近いことがわかる。

Ridge回帰

sklearn.linear_model.Ridge
スクリーンショット 2021-07-14 10.50.48.png

決定係数が 0.819 程度となり、グラフでも真の関数に近いことがわかる。

多項式回帰

sklearn.preprocessing.PolynomialFeatures
スクリーンショット 2021-07-14 10.50.59.png
4~7 次の多項式で決定係数が 0.837 くらいで、8 次以上では 0.840 くらいでほぼ頭打ちの状況となっており、グラフを見ても 4 次以上はすべて真の関数とほぼ重なっていることがわかる。

ラッソ回帰

sklearn.linear_model.Lasso
スクリーンショット 2021-07-14 10.51.12.png
決定係数が 0.821 程度となり、グラフでも真の関数に近いことがわかる。

SVM

sklearn.svm.SVR
スクリーンショット 2021-07-14 10.51.23.png
決定係数が 0.836 程度であり、グラフでも真の関数に近いことがわかる。

NonLiner Regressions via DL by ReLU

from sklearn.model_selection import train_test_split
x_train, x_test, y_train, y_test = train_test_split(data, target, test_size=0.1, random_state=0)
from keras.callbacks import EarlyStopping, TensorBoard, ModelCheckpoint

cb_cp = ModelCheckpoint('/content/drive/My Drive/JDLA_E/Stage2/study_ai_ml_google/skl_ml/out/checkpoints/weights.{epoch:02d}-{val_loss:.2f}.hdf5', verbose=1, save_weights_only=True)
cb_tf  = TensorBoard(log_dir='/content/drive/My Drive/JDLA_E/Stage2/study_ai_ml_google/skl_ml/out/tensorBoard', histogram_freq=0)
def relu_reg_model():
    model = Sequential()
    model.add(Dense(10, input_dim=1, activation='relu'))
    model.add(Dense(1000, activation='relu'))
    model.add(Dense(1000, activation='relu'))
    model.add(Dense(1000, activation='relu'))
    model.add(Dense(1000, activation='relu'))
    model.add(Dense(1000, activation='relu'))
    model.add(Dense(1000, activation='relu'))
    model.add(Dense(1000, activation='relu'))
    model.add(Dense(1000, activation='linear'))
    model.add(Dense(1))

    model.compile(loss='mean_squared_error', optimizer='adam')
    return model
from keras.models import Sequential
from keras.layers import Input, Dense, Dropout, BatchNormalization
from keras.wrappers.scikit_learn import KerasRegressor

# use data split and fit to run the model
estimator = KerasRegressor(build_fn=relu_reg_model, epochs=100, batch_size=5, verbose=1)

history = estimator.fit(x_train, y_train, callbacks=[cb_cp, cb_tf], validation_data=(x_
test, y_test))
Epoch 1/100
18/18 [==============================] - 3s 112ms/step - loss: 1.5787 - val_loss: 1.6068

Epoch 00001: saving model to /content/drive/My Drive/JDLA_E/Stage2/study_ai_ml_google/skl_ml/out/checkpoints/weights.01-1.61.hdf5
Epoch 2/100
18/18 [==============================] - 1s 75ms/step - loss: 1.4145 - val_loss: 1.1353

Epoch 00002: saving model to /content/drive/My Drive/JDLA_E/Stage2/study_ai_ml_google/skl_ml/out/checkpoints/weights.02-1.14.hdf5
Epoch 3/100
18/18 [==============================] - 1s 74ms/step - loss: 1.1486 - val_loss: 1.0478

Epoch 00003: saving model to /content/drive/My Drive/JDLA_E/Stage2/study_ai_ml_google/skl_ml/out/checkpoints/weights.03-1.05.hdf5
Epoch 4/100
18/18 [==============================] - 1s 77ms/step - loss: 0.7643 - val_loss: 0.8383

   ((( 途中の出力は省略 )))

Epoch 00096: saving model to /content/drive/My Drive/JDLA_E/Stage2/study_ai_ml_google/skl_ml/out/checkpoints/weights.96-0.40.hdf5
Epoch 97/100
18/18 [==============================] - 2s 95ms/step - loss: 0.4477 - val_loss: 0.1649

Epoch 00097: saving model to /content/drive/My Drive/JDLA_E/Stage2/study_ai_ml_google/skl_ml/out/checkpoints/weights.97-0.16.hdf5
Epoch 98/100
18/18 [==============================] - 1s 81ms/step - loss: 0.3789 - val_loss: 0.2041

Epoch 00098: saving model to /content/drive/My Drive/JDLA_E/Stage2/study_ai_ml_google/skl_ml/out/checkpoints/weights.98-0.20.hdf5
Epoch 99/100
18/18 [==============================] - 1s 80ms/step - loss: 0.3799 - val_loss: 0.2048

Epoch 00099: saving model to /content/drive/My Drive/JDLA_E/Stage2/study_ai_ml_google/skl_ml/out/checkpoints/weights.99-0.20.hdf5
Epoch 100/100
18/18 [==============================] - 2s 91ms/step - loss: 0.2293 - val_loss: 0.2071

Epoch 00100: saving model to /content/drive/My Drive/JDLA_E/Stage2/study_ai_ml_google/skl_ml/out/checkpoints/weights.100-0.21.hdf5

スクリーンショット 2021-07-14 10.51.37.png
グラフを確認すると真の関数に近いことがわかる。

Lasso 回帰の回帰係数

以下の通り 100 個の回帰係数が表示された。(説明変数が 100 個だから)
スクリーンショット 2021-07-14 10.53.02.png

3.ロジスティック回帰モデル

分類問題(クラス分類)

ある入力(数値)からクラスに分類する問題

識別的アプローチ
$P(C_k \mid \boldsymbol{x})$ を直接モデル化する
識別関数の構成もある(SVMなど)

\left\{
\begin{array}{lll}
f(\boldsymbol{x})>0 & \Rightarrow & C=1\\
f(\boldsymbol{x})<0 & \Rightarrow & C=0
\end{array}
\right.

※ ロジスティック回帰は識別的アプローチ

生成的アプローチ
$P(C_k)$ と $P(\boldsymbol{x} \mid C_k)$ をモデル化し、ベイズの定理を用いて $P(C_k \mid \boldsymbol{x})$ を求める

※ 生成的アプローチを使うと、外れ値の検出や新たなデータを生成(GANなど)することが可能となる。

\begin{align}
P(C_k \mid \boldsymbol{x}) &= \dfrac{P(C_k, \boldsymbol{x})}{P(\boldsymbol{x})} \\
&= \dfrac{P(\boldsymbol{x} \mid C_k)~P(C_k)}{P(\boldsymbol{x})}
\end{align}

シグモイド関数

\sigma(x)=\dfrac{1}{1+e^{-x}}=\dfrac{e^{x}}{1+e^{x}}~
\sigma(x)=\dfrac{1}{1+\exp(-ax)}
  • 入力は実数、出力は必ず 0 〜 1 の値
  • クラス 1 に分類される確率を表現
  • パラメータが変わるとシグモイド関数の形が変わる
    • a を増加させると、x=0 付近での曲線の勾配が増加
    • a を極めて大きくすると、単位ステップ関数( x<0 で f(x)=0、x>0 で f(x)=1 となるような関数)に近づきます
    • バイアス変化は段差の位置
a=2 a=0.5 a=1 a=10
単位ステップ関数に近づいている

シグモイド関数の性質

  • シグモイド関数の微分は、シグモイド関数自身で表現することが可能

連鎖率を使って計算すると・・・

\begin{align}
\dfrac{\partial\sigma(x)}{\partial{x}} 
&= \dfrac{\partial}{\partial{x}} \bigg( \dfrac{1}{1+\exp(-ax)} \bigg) \\
&= \dfrac{\partial}{\partial{x}} \{ 1+\exp(-ax) \}^{-1} \\
&= (-1) \cdot \{ 1+\exp(-ax) \}^{-2} \cdot \exp(-ax) \cdot (-a) \\
&= \dfrac{a\exp(-ax)}{\{ 1+\exp(-ax) \}^2} \\
&= \dfrac{a}{1+\exp(-ax)} \cdot \dfrac{1+\exp(-ax)-1}{1+\exp(-ax)} \\
&= a\sigma(x)(1-\sigma(x))
\end{align}
  • 尤度関数の微分を行う際にこの事実を利用すると計算が容易

モデル数式と幾何学的意味

\large \underset{Y=1となる確率~~~~~~~}{\underline{P(Y=1 \mid \boldsymbol{x})}} 
 = \large \sigma(w_0+w_1x_1+\cdots+w_mx_m)
 = \large \sigma(\boldsymbol{w}^T\boldsymbol{x})
予測:$0$ 予測:$1$
$P(Y=1 \mid \boldsymbol{x})<0.5$ $P(Y=1 \mid \boldsymbol{x})\geq0.5$

最尤推定

ロジスティック回帰モデルの最尤推定

  • 重みパラメータ $\boldsymbol{w}$ を推定
P(Y=y_1 \mid \boldsymbol{x}_1) = p_1^{y_1} (1-p_1)^{1-y_1}
 = \sigma(\boldsymbol{w}^T\boldsymbol{x}_1)^{y_1}
 ( 1 - \sigma(\boldsymbol{w}^T\boldsymbol{x}_1) )^{1-y_1} \\

P(Y=y_2 \mid \boldsymbol{x}_2) = p_2^{y_2} (1-p_2)^{1-y_2}
 = \sigma(\boldsymbol{w}^T\boldsymbol{x}_2)^{y_2}
 ( 1 - \sigma(\boldsymbol{w}^T\boldsymbol{x}_2) )^{1-y_2} \\

\vdots \\

P(Y=y_n \mid \boldsymbol{x}_n) = p_n^{y_n} (1-p_n)^{1-y_n}
 = \sigma(\boldsymbol{w}^T\boldsymbol{x}_n)^{y_n}
 ( 1 - \sigma(\boldsymbol{w}^T\boldsymbol{x}_n) )^{1-y_n} \\

尤度関数

  • $y_1$ 〜 $y_n$ のデータが得られた際の尤度関数
  • 確率変数が独立を仮定▶確率の積に分解可能
  • 尤度関数 $L(\boldsymbol{w})$ はパラメータ $\boldsymbol{w}$ のみに依存する関数
\begin{align}
L(\boldsymbol{w})
 &= \prod_{i=1}^n P(Y=y_i \mid \boldsymbol{x}_i) \\
 &= \prod_{i=1}^n p_i^{y_i} (1-p_i)^{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} \\
\end{align}

交差エントロピー誤差関数

 「尤度関数 $L(\boldsymbol{w})$ を最大化するパラメータを推定(最尤推定)」
$\Longleftrightarrow$「対数尤度関数 $\log L(\boldsymbol{w})$ を最大化するパラメータを推定」
   ※ 同時確率の積が和に、指数が積に変換可能となり、微分が楽になる
$\Longleftrightarrow$「対数尤度関数にマイナスをかけた $-\log L(\boldsymbol{w})$ を最小化するパラメータを推定」
   ※ 誤差関数として使うので、最小二乗法のように最小化するようにする

こうして考えられた $-\log L(\boldsymbol{w})$ が、__交差エントロピー誤差関数 $E(\boldsymbol{w})$__である

\begin{align}
E(\boldsymbol{w})
 &= -\log L(\boldsymbol{w}) \\
 &= -\sum_{i=1}^n \{ y_i\log p_i + (1-y_i)\log(1-p_i) \}
\end{align}

この交差エントロピー誤差関数を最小化することが、もともとの尤度関数を最大化(最尤推定)することになる

勾配降下法

  • 反復学習によりパラメータを逐次的に更新するアプローチの一つ
  • ロジスティック回帰モデルでの最尤法に関して、対数尤度関数をパラメータで微分して $0$ になる値を解析的に求めることは困難であるため勾配降下法を使う
    • 線形回帰モデルの MSE では、パラメータに関する微分が0になる値を解析的に求めることが可能
  • $\eta$ は学習率と呼ばれるハイパーパラメータでモデルのパラメータの収束しやすさを調整
\begin{align}
\boldsymbol{w}^{(k+1)}
 &= \boldsymbol{w}^{(k)}
  - \eta\dfrac{\partial E(\boldsymbol{w})}{\partial \boldsymbol{w}} \\
 &= \boldsymbol{w}^{(k)}
  + \eta \sum_{i=1}^{n} (y_i-p_i)\boldsymbol{x}_i \\
\end{align}

__「パラメータが更新されなくなった」ということは、「勾配が $0$ になった」__ということ
(少なくとも反復学習で探索した範囲では最適な解がもとめられたことになる。)

勾配降下法の問題点と解決策

  • パラメータを更新するのに $n$ 個全てのデータに対する和を求める必要がある
    (計算リソース)nが巨大になった時にデータをオンメモリに載せる容量が足りない
    (計算コスト)計算時間が莫大になる

  • __確率的勾配降下法__を利用して解決

確率的勾配降下法(SGD)

  • データを一つずつランダムに(「確率的」に)選んでパラメータを更新
  • 勾配降下法でパラメータを $1$ 回更新するのと同じ計算量でパラメータを $n$ 回更新できるので、効率よく最適な解を探索可能
\boldsymbol{w}^{(k+1)}
 = \boldsymbol{w}^{(k)}
  + \eta (y_i-p_i)\boldsymbol{x}_i \\
勾配降下法 確率的勾配降下法(SGD)
$\boldsymbol{w}^{(k+1)} = \boldsymbol{w}^{(k)} + \eta \sum_{i=1}^{n} (y_i-p_i)\boldsymbol{x}_i$ $\boldsymbol{w}^{(k+1)} = \boldsymbol{w}^{(k)} + \eta (y_i-p_i)\boldsymbol{x}_i$

モデルの評価

混同行列(confusion matrix)

について説明○各検証データに対するモデルの予測結果を4つの観点(表)で分類し、それぞれに当てはまる予測結果の個数をまとめた表
スクリーンショット 2021-07-14 22.09.30.png

再現率(Recall)と適合率(Precision)

再現率(Recall)

\mathrm{Recall}=\frac{TP}{TP+FN}
  • 「本当にPositiveなもの」の中からPositiveと予測できる割合
    • NegativeなものをPositiveとしてしまう事象については考えていない
  • 「誤り(False Positive)が多少多くても抜け漏れは少ない」予測をしたい際に利用
    • 使用例)病気の検診で「陽性であるものを陰性と誤診(False Negative)」としてしまうのを避けたい。「陰性を陽性であると誤診(False Positive)」とするものが少し増えたとしても再検査すればいい。

適合率(Precision)

\mathrm{Precision}=\frac{TP}{TP+FP}
  • モデルが「Positiveと予測」したものの中で本当にPositiveである割合
    • 本当にPositiveなものをNegativeとしてしまう事象については考えていない
  • 見逃し(False Negative)が多くてもより正確な」予測をしたい際に利用
    • 使用例)「重要なメールをスパムメールと誤判別」されるより「スパムと予測したものが確実にスパム」である方が便利。スパムメールを検出できなくても(False Negative)、自分でやればいい。

F値

\mathrm{F}値 = \frac{2\mathrm{Recall}\cdot\mathrm{Precision}}{\mathrm{Recall}+\mathrm{Precision}}
  • 理想的にはどちらも高いモデルがいいモデルだが、両者はトレードオフの関係にあり、どちらかを小さくすると、もう片方の値が大きくなってしまう。
  • Precision と Recall の調和平均
    • Recall と Precision のバランスを示している。
    • 高ければ高いほど Recall と Precision がともに高くなる。

<参考動画>

実装演習(タイタニックデータ解析)

0. データ表示

  • 各種モジュールをインポートする。

  • タイタニックのCSVファイルを読み込む。

  • データセットの先頭5つのデータを確認する。

  • すべてのデータを表示して、データ数と属性数を確認する。

    出力の最後の行をみると、データ数が891で属性数が12あることがわかる。

1. ロジスティック回帰

1-1. 不要なデータの削除・欠損値の補完

  • 予測に必要ないと思われる説明変数を取り除く。

    ここでは、4つの説明変数(乗客ID、乗客名、チケットNo.、船室)をドロップしているが、船室は船の中の位置(船首側が先に沈んでしまい船室から逃げる時間がなかったなど)などの理由で生死に関係しているかもしれないとか考えたりする・・・。

  • 欠損値を含むデータを表示する。

    ここでは、年齢に欠損のあるデータが存在することがわかった。

  • 年齢Ageに欠損があるデータに対して、平均の値で補完し、それを新たな説明変数AgeFillに入れる。

    ここで、説明変数Ageと新しくつくった説明変数Agefillをデータの先頭30個で比較してみた。

Age AgeFill

もともとAgeのデータがあったデータのAgeFillにはもとのAgeの値が、Ageに欠損のあったデータのAgeFillにはAgeの平均値の29.699118が入っていることが確認できた。

1-2. 実装(チケット価格から生死を判別)

  • 運賃だけのリストをdata1として作成する。(運賃だけが説明変数)

  • 生死フラグだけのリストをlabel1として作成する。(生死フラグが目的変数)

  • ロジスティック回帰モデルのモジュールをインポートする。

  • ロジスティック回帰モデルのインスタンスとしてmodelを作成する。

  • ロジスティック回帰モデルで重みを学習する。

  • 説明変数data1(運賃Fare)の値を62として、生き残ったかどうか確認する。

    1が出力されたので生き残ったことがわかった。

  • 説明変数data1(運賃Fare)の値を62として、各クラスに属する確率を予測する。

    ここでは死亡(0)の確率が0.49978123、生存(1)の確率が0.50021877となっており、運賃が62払っている人の生存確率がほぼ半々という予測がでている。
    ここで、運賃62の半分の運賃31の場合と2倍の運賃124の場合でも予測してみた。

運賃 31 運賃 62 運賃 124
死ぬ確率が約6割 生きるか死ぬか半々 生き残れる確率約7割

高い運賃を払っていた人の方が生存確率が高かったという予測結果となった。

  • 切片と重みを確認する。

    ここでは、切片(intercept_)が -0.94131796、重み(coef_)が 0.01519666 となっていることが確認できた。今回は説明変数がひとつなので重みもひとつ表示されている。

  • 学習したモデルに説明変数を入力し、予測値として出力される生死の確率をグラフで描画する。

    オレンジのドットが死亡(0)する確率、緑のドットが生存(1)する確率、赤い実線はシグモイド関数である。緑の生存する確率がシグモイド関数に一致していることがわかる。

1-3. 実装(2変数から生死を判別)

  • 新しく説明変数 Gender を作成し、説明変数 Sex が female であれば 0、male であれば 1 を設定する。
    スクリーンショット 2021-07-16 15.00.00.png

  • 説明変数 Gender が追加され、値が入っていることを確認する。
    スクリーンショット 2021-07-16 15.00.09.png

  • もうひとつ新しい説明変数 Pclass_Gender を作成し、乗客の社会階級を示す Pclass と 先ほど作った Gender の合計値を設定する。
    スクリーンショット 2021-07-16 15.00.20.png
    これが意味するのは、Pclass_Gender は値が小さいほど社会階級が高いということと、Gender はレディーファーストなどの文化で女性が優先されるなどの傾向があると思われ女性の値の方を小さい値にしたということを利用して、Pclass と Gender の合計である Pclass_Gender も小さいほど優先されるという仮定に基づく設定となっている。

  • 説明変数 Pclass_Gender が追加され、値が入っていることを確認する。
    スクリーンショット 2021-07-16 15.00.30.png

  • 不要となった説明変数をドロップする。
    スクリーンショット 2021-07-16 15.00.39.png

  • 最終的に残った説明変数を確認する。
    スクリーンショット 2021-07-16 15.00.49.png

  • 横軸が AgeFill、縦軸が Pclass_Gender として、生死のどちらだったかをグラフに描画する。(生存が青、死亡が赤である。)
    スクリーンショット 2021-07-16 15.01.15.png
    Pclass_Gender が 1 ということは、階級がもっとも高い女性ということであり、ほぼ生存していることがわかる。また Pclass_Gender が 4 となる階級が最も低い男性はかなりの割合で死亡していることがわかる。

  • AgeFill、Pclass_Gender の2つだけのリストを data2 として作成する。(説明変数)
    スクリーンショット 2021-07-16 15.01.33.png

  • data2 が作られていることを確認する。
    スクリーンショット 2021-07-16 15.01.43.png

  • 生死フラグだけのリストを label2 として作成する。(目的変数)
    スクリーンショット 2021-07-16 15.01.53.png

  • ロジスティック回帰モデルのインスタンスとして model2 を作成する。
    スクリーンショット 2021-07-16 15.02.03.png

  • ロジスティック回帰モデルで重みを学習する。
    スクリーンショット 2021-07-16 15.02.12.png

  • 説明変数data2の要素のAgeFillを10、Pclass_Genderを1として、生き残ったかどうか確認する。
    スクリーンショット 2021-07-16 15.02.23.png
    1が出力されたので生き残ったことがわかった。

  • 説明変数data2の要素のAgeFillを10、Pclass_Genderを1として、各クラスに属する確率を予測する。
    スクリーンショット 2021-07-16 15.02.33.png
    ここでは死亡(0)の確率が0.03754749、生存(1)の確率が0.96245251となっており、階級が高い10歳の女の子の生存確率が0.96245251と非常に高いことがわかった。

<境界線の式>
$~w_1x + w_2y + w_0 = 0~\Longleftrightarrow~y = \dfrac{-w_1x - w_0}{w_2}~$

  • 上記の境界線の式を使って、グラフに境界線をプロットする。
  • $~x:AgeFill$、$y:Pclass\_Gender$、$w_0:切片$、$w_1,w_2:重み~$
    スクリーンショット 2021-07-16 15.02.54.png

2. モデル評価

ここは実行結果のキャプチャのみ

2-1. 混同行列とクロスバリデーション

スクリーンショット 2021-07-16 15.03.11.png
スクリーンショット 2021-07-16 15.03.20.png
スクリーンショット 2021-07-16 15.03.29.png
スクリーンショット 2021-07-16 15.03.39.png
スクリーンショット 2021-07-16 15.03.49.png
スクリーンショット 2021-07-16 15.03.57.png
スクリーンショット 2021-07-16 15.04.07.png
スクリーンショット 2021-07-16 15.04.16.png
スクリーンショット 2021-07-16 15.04.24.png
スクリーンショット 2021-07-16 15.04.34.png
スクリーンショット 2021-07-16 15.04.52.png
スクリーンショット 2021-07-16 15.05.04.png
スクリーンショット 2021-07-16 15.05.12.png
スクリーンショット 2021-07-16 15.05.21.png
スクリーンショット 2021-07-16 15.05.36.png
スクリーンショット 2021-07-16 15.05.50.png
スクリーンショット 2021-07-16 15.06.06.png
スクリーンショット 2021-07-16 15.06.17.png

4.主成分分析(PCA)

  • 多変量データの持つ構造をより少数個の指標に圧縮
    • 変量の個数を減らすことに伴う、情報の損失はなるべく小さくしたい
    • 少数変数を利用した分析や可視化(2・3次元の場合)が実現可能

<参考動画>

学習データ
$m$ 次元の説明変数

\boldsymbol{x}_i=(x_{i1}, x_{i2}, \cdots, x_{im}) \in\mathbb{R}^m

平均(ベクトル)
$n$ 個の学習データの平均

\bar{\boldsymbol{x}}=\frac{1}{n} \sum_{i=1}^n \boldsymbol{x}_i

データ行列
中心化した学習データを並べてつくった$n\times m$行列(原点の周りに分散するようにしたということ)

\bar{X}
 =(\boldsymbol{x}_1-\bar{\boldsymbol{x}}, \cdots, \boldsymbol{x}_n-\bar{\boldsymbol{x}})^T
 =
 \begin{pmatrix}
  \boldsymbol{x}_1-\bar{\boldsymbol{x}} \\
  \vdots \\
  \boldsymbol{x}_n-\bar{\boldsymbol{x}}
 \end{pmatrix} \in\mathbb{R}^{n\times m}

分散共分散行列

\Sigma=Var(\bar{X})=\frac{1}{n}\bar{X}^T\bar{X}

線形変換後のベクトル
$~j~$は射影軸のインデックス(圧縮先の次元が$1$次元なら $j-1$、$2$次元なら $j=2$など)

\underset{n\times 1}{\boldsymbol{s}_j}
 =(s_{1j}, \cdots, s_{nj})^T
 =\begin{pmatrix}
  s_{1j} \\
  \vdots \\
  s_{nj}
 \end{pmatrix}
 =\underset{n\times m}{\bar{X}}\underset{m\times 1}{\boldsymbol{a}_j}\qquad \boldsymbol{a}_j
 \in\mathbb{R}^m \\
  • 係数ベクトルが変われば線形変換後の値が変化
    • 情報の量を分散の大きさと捉える
    • 線形変換後の変数の分散が最大となる射影軸を探索

線形変換後の分散

~\begin{align}
Var(\boldsymbol{s}_j)
 &= \frac{1}{n} \boldsymbol{s}_j^T\boldsymbol{s}_j
  = \frac{1}{n} (\bar{X}\boldsymbol{a}_j)^T(\bar{X}\boldsymbol{a}_j)
  = \frac{1}{n} \boldsymbol{a}_j^T~\bar{X}^T \bar{X}~\boldsymbol{a}_j \\
 &= \boldsymbol{a}_j^T Var(\bar{X})~\boldsymbol{a}_j
\end{align}~

ラグランジュの未定乗数法

  • 以下の制約付き最適化問題を解く
    • ノルムが1となる制約を入れる(制約を入れないと無限に解がある)
\overset{目的関数}{\arg \underset{\boldsymbol{a}\in\mathbb{R}^m}{\max} {\boldsymbol{a}_j}^T Var(\bar{X}) \boldsymbol{a}_j} \qquad
\overset{制約条件}{{\boldsymbol{a}_j}^T \boldsymbol{a}_j=1}
  • 制約つき最適化問題の解き方
    • ラグランジュ関数$~E(\boldsymbol{a}_j)~$を最大にする係数ベクトルを探索(微分して0になる点)
\underset{ラグランジュ関数}{E(\boldsymbol{a}_j)}
 = \underset{目的関数}{\underline{{\boldsymbol{a}_j}^T Var(\bar{X}) \boldsymbol{a}_j}}
 - \overset{ラグランジュ乗数}{\lambda} \underset{制約条件}{\underline{{(\boldsymbol{a}_j}^T \boldsymbol{a}_j-1)}}
  • ラグランジュ関数を微分して最適解を求める
\frac{\partial E(\boldsymbol{a}_j)}{\partial\,\boldsymbol{a}_j}
= 2Var(\bar{X})~\boldsymbol{a}_j-2\lambda~\boldsymbol{a}_j=0 \\
\therefore Var(\bar{X})~\boldsymbol{a}_j = \lambda~\boldsymbol{a}_j

  この解は、固有値と固有ベクトルの定義そのものである
  つまり、元のデータの分散共分散行列$~Var(\bar{X})~$の固有値$~\lambda~$と固有ベクトル$~\boldsymbol{a}_j~$が、上記の制約付き最適化問題の解となる

  射影先の分散は固有値と一致

\begin{align}
Var(\boldsymbol{s}_1)
 &= \boldsymbol{a}_1^T~Var(\bar{X})~~\boldsymbol{a}_1 \\
 &= \lambda_1~\boldsymbol{a}_1^T~\boldsymbol{a}_1 \\
 &= \lambda_1
\end{align}

<参考動画>

  • 分散共分散行列は半正定値対称行列 ▶ 固有値は必ず 0 以上・固有ベクトルは直交

k 番目の固有値を昇順に並べ、対応する固有ベクトルを第 k 主成分と呼ぶ。

寄与率

  • 第1~元次元分の主成分の分散は、元のデータの分散と一致
    • 2次元のデータを2次元の主成分で表示した時、固有値の和と元のデータの分散が一致
    • 第k主成分の分散は主成分に対応する固有値
\underset{元データの総分散}{V_{total}}
 = \underset{主成分の総分散}{\sum_{i=1}^m \lambda_i}
  • 寄与率:第k主成分の分散の全分散に対する割合(第k主成分が持つ情報量の割合)
\underset{寄与率}{c_k}
= \frac{\overset{第k主成分の分散}{\lambda_k}}{~~\underset{主成分の総分散}{\displaystyle\sum_{i=1}^m \lambda_i}~~}
  • 累積寄与率:第1-k主成分まで圧縮した際の情報損失量の割合
\underset{累積寄与率}{r_k}
= \frac{\overset{第1〜k主成分の分散}{\displaystyle\sum_{j=1}^k \lambda_j}}
       {~~\underset{主成分の総分散}{\displaystyle\sum_{i=1}^m \lambda_i}~~}

実装演習(乳がんデータ分析)

  • 各種モジュールをインポートする。
    スクリーンショット 2021-07-16 18.57.06.png

  • 乳がんのCSVファイルを読み込む。
    スクリーンショット 2021-07-16 18.57.42.png

  • 乳がんデータの行数と列数を確認する。
    スクリーンショット 2021-07-16 18.58.12.png
    このデータは569行、33列であることがわかる。

  • データを表示し確認する。
    スクリーンショット 2021-07-16 18.58.29.png
    横にスクロールして確認すると、最後の33列目に Unnamed: 32 というNaNで埋め尽くされたようなゴミデータが見つかる。
    スクリーンショット 2021-07-16 18.58.49.png

  • Unnamed: 32 をドロップし、再度表示して確認する。
    スクリーンショット 2021-07-16 18.59.01.png
    569 rows × 32 columns と表示され、1列減っていることが確認できる。
    スクリーンショット 2021-07-16 18.59.28.png
    また右に横スクロールして、最後の Unnamed: 32 の列がなくなっていることも確認できる。

説明変数は3列以降、目的変数は2列目のdiagnosisとしロジスティック回帰で分類(diagnosis: 診断結果(良性がB/悪性がM))

  • 目的変数 diagnosis の値が M ならば 1、それ以外なら 0 を y に入れる。
    スクリーンショット 2021-07-16 18.59.46.png

  • 説明変数の radius_mean 以降の全ての列を X に入れる。(id, diagnosisは含まれないので 30 次元)
    スクリーンショット 2021-07-16 18.59.58.png

  • ロジスティック回帰で学習を行う。
    スクリーンショット 2021-07-16 19.01.02.png
    結果はでているようであるが、同じWarningが繰り返し表示される。
    どうやら「イテレーションの限界に達したけど、収束しないから max_iter の値を増やせ」とのことらしかったので、10行目の LogisticRegressionCV の最後に max_iter=200(デフォルトは 100 らしい)をつけてやりなおしてみた。
    スクリーンショット 2021-07-16 19.01.36.png
    すると、無事 Warning が消えて結果のみが出力されるようになった。
    訓練スコアが 0.988、検証スコアが 0.972。

  • 次元圧縮する前の 30 次元で寄与率をグラフ化し確認する。(縦軸が寄与率、横軸が第 n 主成分)
    スクリーンショット 2021-07-16 19.02.06.png
    ざっと見ると、第1主成分(40%強)、第2主成分(20%弱)、第3主成分(10%弱)、第4主成分(5%強)なので、第1〜4主成分の寄与率の合計が80%くらいとなっていることがわかる。

  • 2 次元まで次元圧縮してプロットする。
    スクリーンショット 2021-07-16 19.02.16.png
    重なってしまっているところはあるが、いい感じで分離できているのではないだろうか。

5.アルゴリズム

k 近傍法(kNN)

  • 分類問題のための機械学習手法

k 近傍法(kNN)のアルゴリズム

  • 最近傍のデータを k 個取ってきて、それらが最も多く所属するクラスに識別

k の変化による結果への影響

  • k を変化させると結果も変わる
k=1(最近傍法) k=3 k=5
紫(1), 黄(0)
→ 紫クラスに分類
紫(2),黄(1)
→ 紫クラスに分類
紫(2),黄(3)
→ 黄クラスに分類
  • k を大きくすると決定境界は滑らかになる
k=1(最近傍法) k=3 k=10

<参考動画>

実装演習

  • 各種モジュールをインポートする。
    スクリーンショット 2021-07-16 23.13.54.png

訓練データ生成

  • 訓練データを生成する。
    スクリーンショット 2021-07-16 23.14.20.png

  • 生成した訓練データをプロットして確認する。
    スクリーンショット 2021-07-16 23.14.32.png

学習

  • 陽に訓練ステップはない。

予測

  • numpy 実装
    スクリーンショット 2021-07-16 23.14.54.png
    スクリーンショット 2021-07-16 23.15.03.png

  • scikit-learn 実装
    スクリーンショット 2021-07-16 23.15.14.png
    スクリーンショット 2021-07-16 23.36.52.png

k 平均法(k-means)

  • 教師なし学習
  • クラスタリング手法(特徴の似ているもの同士をグループ化)
  • 与えられたデータを k 個のクラスタに分類する

k-means(k 平均法)のアルゴリズム

k-meansアルゴリズムは以下のとおりである。

  1. 各クラスタ中心の初期値を設定する
  2. 各データ点に対して各クラスタ中心との距離を計算し、最も距離が近いクラスタを割り当てる
  3. 各クラスタの平均ベクトル(中心)を計算する
  4. 収束するまで2, 3の処理を繰り返す

手順①

各クラスタ中心の初期値を設定する。

手順②

各データ点に対して、各クラスタ中心との距離を計算し、最も距離が近いクラスタを割り当てる。

手順③

各クラスタの平均ベクトル(中心)を計算する。

手順④

クラスタの再割り当てと、中心の更新を繰り返す。

初期値や k の変化による結果への影響

  • 中心の初期値を変えるとクラスタリング結果も変わりうる
    ※ 初期値の問題を改善する、k-means++ というものがある。

  • k の値を変えるとクラスタリング結果も変わる

<参考動画>

実装演習

スクリーンショット 2021-07-16 23.38.49.png

データ生成

スクリーンショット 2021-07-16 23.39.53.png

学習

スクリーンショット 2021-07-16 23.40.40.png

クラスタリング結果

  • numpy 実装
    スクリーンショット 2021-07-16 23.40.50.png

  • scikit-learn 実装
    スクリーンショット 2021-07-16 23.41.20.png

6.サポートベクターマシーン(SVM)

汎化性能が高さ、応用分野の広さ(回帰問題や教師なし学習など)からデータ分析で近年もっとも注目を集めているモデル。

2クラス分類

与えられた入力データが2つのカテゴリー(クラス)のどちらに属するかを識別する問題を__2クラス分類問題__という。言い換えると、$~\boldsymbol{x}~$(特徴ベクトル、入力ベクトル)を与えたときに予測値$~y~$(ラベル)を返す関数をつくる問題である。

分類問題のためのSVMを他の問題用と区別して__サポートベクトル(SV)分類__ と呼ぶ。

決定関数と分類境界

特徴ベクトル$~\boldsymbol{x}~$がどちらのクラスに属するかを判定するときに使われる関数$~f(\boldsymbol{x})~$を__決定関数__という。

f(\boldsymbol{x}) = \boldsymbol{w}^T\boldsymbol{x} + b
  • $~\boldsymbol{w}~$は特徴ベクトル$~\boldsymbol{x}~$と同じ次元の数値ベクトル

ある入力データ$~\boldsymbol{x}~$に対して決定関数$~f(\boldsymbol{x})~$を計算し、その符号により2つのクラスに分類する。

y = \mathrm{sgn}~f(\boldsymbol{x}) =
\left\{
\begin{array}{ll}
+1 & (f(\boldsymbol{x}) > 0) \\
-1 & (f(\boldsymbol{x}) < 0)
\end{array}
\right.
  • $~y~$はラベル
  • $~\mathrm{sgn}(\cdot)~$は符号関数:引数が正の場合は+1、負の場合には−1を返す・

具体的なイメージを掴むため、特徴ベクトルが2次元の場合を考える。
特徴ベクトルが2次元のときの決定関数は、

\begin{align}
f(\boldsymbol{x})
 &= \boldsymbol{w}^T\boldsymbol{x} + b \\
 &= (w_1, w_2)
 \begin{pmatrix}
  x_1 \\
  x_2
 \end{pmatrix} + b \\
 &= w_1x_1 + w_2x_2 + b

\end{align}

となり、決定関数$~f(\boldsymbol{x})~$は下図に示すように平面を表す方程式となる。
(点は特徴ベクトル$~\boldsymbol{x} = (x_1, x_2)^T~$)

スクリーンショット 2021-07-17 12.07.20.png

この図で、$~f(\boldsymbol{x})>0~$の場合(青の領域) では$~y=1~$、$~f(\boldsymbol{x})<0~$の場合(赤の領域) では$~y=−1~$ となる。

$~x_1-x_2~$平面と$~f(\boldsymbol{x})~$の平面の交線$~f(\boldsymbol{x})=0~$が、特徴ベクトル$~x=(x_1, x_2)^T~$を2つのクラスに分ける境界線になっていることが分かる。一般にこのような境界を__分類境界__という。(2次元データの場合には分類境界は「直線」、n 次元データの場合には「超平面」)

すなわち、決定関数$~f(\boldsymbol{x})~$のパラメータ$~\boldsymbol{w}, b~$を決定することが、適切な分類境界を決定することとなる。

線形サポートベクトル分類

ハードマージン

訓練データは特徴ベクトルとラベルのセットである。

(\boldsymbol{x}_i, y_i) \quad (i = 1, 2, · · · , n)

n 個の訓練データをすべて正しく分類できる$~\boldsymbol{w}~$と$~b~$の組が存在する場合(__分離可能__という)、決定関数によるラベルの予測値 $~\mathrm{sgn}~f(\boldsymbol{x})~$と訓練データの値$~y_i~$の符号が一致しているということなので次式が成り立つ。

y_i~f(\boldsymbol{x}_i)>0 \quad (i = 1, 2, · · · , n)

もし正しく分類できていないデータがある場合には、そのデータ$~(\boldsymbol{x}_i, y_i)~$に対して$~\mathrm{sgn}~f(\boldsymbol{x})~$と$~y_i~$の符号が逆になっているため$~y_i~f(\boldsymbol{x}_i)<0~$となる。

一般に訓練データを分離できる分類境界は複数存在する。SVMでは、それぞれのクラスのデータが分類境界からなるべく離れるようにして「最適な」分類境界を決定する。

分類境界を挟んで2 つのクラスがどのくらい離れているかを__マージン__という。

スクリーンショット 2021-07-17 13.03.25.png

SVM ではマージンが大きいほど良い分類境界と考えるため、なるべく大きなマージンを持つ分類境界を探すことになる。これを__マージン最大化__という。

マージンを最大化することは、分類境界$~f(\boldsymbol{x})=0~$と分類境界から最も近くにあるデータ$~\boldsymbol{x}_i~$との距離を最大化することになる。最終的に SVM の目的関数は次のように書ける。

\underset{\boldsymbol{w}, b}{\max}\bigg[ \underset{i}{\min} \frac{y_i[\boldsymbol{w}^T\boldsymbol{x}_i+b]}{||\boldsymbol{w}||} \bigg]
=
\underset{\boldsymbol{w}, b}{\max} \frac{M(\boldsymbol{w}, b)}{||\boldsymbol{w}||}

この目的関数の構造は、次のように2つの要素からなっている。

\begin{align}
&\underset{i}{\min} \big[ y_i[\boldsymbol{w}^T\boldsymbol{x}_i+b] \big] = M(\boldsymbol{w}, b) \\
&\Leftrightarrow すべての~i~に対して~y_i[\boldsymbol{w}^T\boldsymbol{x}_i+b]\geq M(\boldsymbol{w}, b) \tag{1} \\
\\
&\underset{\boldsymbol{w}, b}{\max}\dfrac{M(\boldsymbol{w}, b)}{||\boldsymbol{w}||} \tag{2}
\end{align}

メインとなる式$~(2)~$の最大化の条件から、決定関数のパラメータ$~\boldsymbol{w}~$と$~b~$が決定される。

$~\frac{M(\boldsymbol{w}, b)}{||\boldsymbol{w}||}~$を最大化したい場合には$~M(\boldsymbol{w}, b)~$は大きい方がよい。しかし$~M(\boldsymbol{w}, b)~$が大きくなりすぎると、式$~(1)~$の条件を満たせなくなる。つまり「式$~(1)~$の条件を満たせる範囲で、式$~(2)~$の最大化を行う」というのが、この目的関数の意味するところである。一般に式$~(1)~$のような条件のことを__制約条件__という。

これを扱いやすい形に変形すると以下となる。

\underset{\boldsymbol{w}, b}{\min} \frac{1}{2} || \boldsymbol{w} ||^2
,\quad ただし、すべての~i~に対して~y_i[\boldsymbol{w}^T\boldsymbol{x}_i+b]\geq 1

この最適化問題は、分離可能性を仮定した場合の線形 SV 分類の標準的な定式化となっている。
なお、分離可能性を仮定した SV 分類のことを一般に__ハードマージン__という。

SV 分類では「分類境界$~f(\boldsymbol{x})=0~$と分類境界から最も近くにあるデータ$~\boldsymbol{x}_i~$との距離の最大化」に基づいて分類境界が決定される。そのため、分類境界に最も近いデータ以外については取り除いてしまっても SV 分類によって得られる分類境界は変化しない。

このように分類境界に最も近いデータ$\boldsymbol{x}_i$が分類境界を支えていると解釈できるため、この$~\boldsymbol{x}_i~$のことを__サポートベクトル__と呼ぶ

ソフトマージン

ハードマージンでは訓練データを完璧に分類できる決定関数$~f(\boldsymbol{x})~$が存在するという仮定をした。しかし現実的にはこの仮定は強すぎるため、この仮定をなくし SV 分類を分離可能でないデータに適用できるようにする方法がもとめられる。この方法を__ソフトマージン__という。

ハードマージンからソフトマージンへの拡張は、ハードマージンの制約条件である「すべての$~i~$に対して$~y_i[\boldsymbol{w}^T\boldsymbol{x}_i+b]\geq 1~$」を多少の分類誤りは許すように緩和することで実現する。

このために新たに非負の変数$~\xi_i \geq 0 (i=1, \cdots, n)~$を導入し、この制約条件を次のように変更する。

y_i[\boldsymbol{w}^T\boldsymbol{x}_i+b]\geq 1-\xi_i \qquad (i=1, \cdots, n)

この$~\xi_i~$はマージン内に入るデータや誤分類されたデータに対する誤差を表す変数となっており、__スラック変数__という。

ハードマージンのときは、分類境界とサポートベクトルとの距離を用いてマージンを定義したが、ソフトマージンでは下図のように$~f(\boldsymbol{x})=1~$と$~f(\boldsymbol{x})=-1~$の間の距離をマージンと解釈する。そして、このマージンを最大化しつつ、分類の誤差$~\xi_i~$を最小化するように分類境界を決定する。

このソフトマージンの場合のSV 分類における最適化問題は、数式的には次のように表現される。

\begin{align}
&\underset{\boldsymbol{w}, b, \boldsymbol{\xi}}{\min} 
\bigg[ \frac{1}{2} || \boldsymbol{w} ||^2 + C\sum_{i=1}^n\xi_i \bigg], \tag{1} \\
&ただし、y_i[\boldsymbol{w}^T\boldsymbol{x}_i+b]\geq 1-\xi_i ~~ (i=1, \cdots, n) \tag{2}
\end{align}

$~\boldsymbol{\xi}=(\xi_1, \cdots, \xi_n)^T~$とする。また、係数$~C~$は__正則化係数__と呼ばれる正の定数で、学習前に値を与えておく必要があるハイパーパラメータである。

目的関数である式$~(1)~$の第1項はハードマージンの場合と同様にマージン最大化の働きをもつ。第2項は本来の制約条件$~y_i[\boldsymbol{w}^T\boldsymbol{x}_i+b]\geq 1~$からのズレを表す$~\xi_i~$がなるべく小さくなるように抑制している。

この項によって、たとえマージンが大きくても誤分類がたくさん発生するような分類境界はつくられにくくなっており、正則化係数$~C~$がこの抑制の度合いを調整するパラメータになっている。

|正則化係数$~C~$|誤分類|マージン|
|:-:|:-:|:-:|:-:|
|小さい
↕︎
大きい|許容されやすくなる
↕︎
少なくなる|ハードマージンから遠のく
↕︎
ハードマージンに近づく|
|$~C \rightarrow \infty~$|誤分類なし|ハードマージンに一致 ($~\xi_i=0~$)|

$~C~$が大きいほど誤分類が少なく良い分類境界を構成できるが、$C~$が大きすぎるとデータが分離可能でない場合、あらゆる分離境界に対して目的関数が発散して計算できなくなる可能性がある。そのため、実装する際にはデータに合わせてちょうど良い$~C~$の値を決定する必要がある。

$~C~$の値の決定には、交差検証法などが用いられる。

SVM における双対表現

線形 SV 分類を用いて分類境界を決定する問題は、以下の最適化問題に帰着した。一般にこれらの最適化問題のことを、SV 分類の__主問題__という。

\begin{align}
ハ&ードマージン \\
&\underset{\boldsymbol{w}, b}{\min} \frac{1}{2} || \boldsymbol{w} ||^2
, \quad y_i[\boldsymbol{w}^T\boldsymbol{x}_i+b]\geq 1 ~~ (i=1, \cdots, n)\\
&\\
ソ&フトマージン \\
&\underset{\boldsymbol{w}, b, \boldsymbol{\xi}}{\min} 
\bigg[ \frac{1}{2} || \boldsymbol{w} ||^2 + C\sum_{i=1}^n\xi_i \bigg], \quad y_i[\boldsymbol{w}^T\boldsymbol{x}_i+b]\geq 1-\xi_i ~~ (i=1, \cdots, n)
\end{align}

主問題を解けば分類境界を決定できるが、以下のメリットがあるため主問題と等価な__双対__(そうつい)__問題__の方を解くことが一般的となっている。

  • 主問題と比べて双対問題の方が変数を少なくできる
  • 分類境界の非線形化を考える上で双対問題の形式(双対形式)の方が有利となる

SVM の双対問題

主問題は$~(\boldsymbol{w}, b, \boldsymbol{\xi})~$というパラメータに関する最適化問題であったが、双対問題は以下のように$~\alpha~$のみに関する最適化問題となっており、変数が減って解くのが簡単になっている。

\begin{align}
ハ&ードマージン \\
&\underset{\alpha}{\max} \bigg[ -\dfrac{1}{2} \sum_{i=1}^n \sum_{j=1}^n \alpha_i \alpha_j y_i y_j \boldsymbol{x}_i^T \boldsymbol{x}_j + \sum_{i=1}^n \alpha_i \bigg], \quad \sum_{i=1}^n \alpha_i y_i = 0, \quad 0\leq\alpha_i ~~ (i=1, \cdots, n) \\
&\\
ソ&フトマージン \\
&\underset{\alpha}{\max} \bigg[ -\dfrac{1}{2} \sum_{i=1}^n \sum_{j=1}^n \alpha_i \alpha_j y_i y_j \boldsymbol{x}_i^T \boldsymbol{x}_j + \sum_{i=1}^n \alpha_i \bigg], \quad \sum_{i=1}^n \alpha_i y_i = 0, \quad 0\leq\alpha_i\leq C ~~ (i=1, \cdots, n) \\
\end{align}

カーネルを用いた非線形分離への拡張

データの分布によっては線形分離ではうまく分類できないケースがある。SVM では特徴ベクトルを非線形変換して、その空間での線形分類を行う「カーネルトリック」と呼ばれる手法により非線形分離を行うことも可能である。

スクリーンショット 2021-07-17 18.16.18.png

上図は線形分離できないデータの例であるが、このままでは単純な直線(超平面) では分類できない。このような場合には、入力データを高次元に拡張することで線形分離が可能になるケースがある。

上図の$~(x_1,~x_2)~$という2次元データに対して、$~(x_1,~x_2,~x_1^2+x_2^2)~$という3次元データをつくると、以下のようにこの3次元空間では平面(超平面)により2つのクラスに分類することが可能になる。

スクリーンショット 2021-07-17 18.16.42.png

このように、入力データをより次元の高い空間(特徴空間)のデータに拡張し、その高次元空間で線形分離を行い、その結果を元の入力データに落とし込むというのが非線形分離の基本的なアイディアである。

このような高次元データへの拡張を数学的には写像(map)という。上の例では次の写像を考えた。

\phi(\boldsymbol{x}) : (x_1,~x_2)~\rightarrow~(x_1,~x_2,~x_1^2+x_2^2)

ここで$~\boldsymbol{x} = (x_1, x_2)^T~$は元の特徴ベクトルを表す。この$~\boldsymbol{x}~$から$~(x_1,~x_2,~x_1^2+x_2^2)~$の高次元データへ写像する関数を$~\phi(\boldsymbol{x})~$としている。これを一般化すれば、入力空間で n 次元であるデータをより高次のr 次元の特徴空間へ変換する関数は次のようになる。

\boldsymbol{\phi}(\boldsymbol{x}) =
\begin{pmatrix}
 \phi_1(\boldsymbol{x}) \\
 \vdots \\
 \phi_r(\boldsymbol{x})
\end{pmatrix}

この写像$~\boldsymbol{\phi}(\boldsymbol{x})~$により、SV 分類における双対問題での目的関数は、次のように特徴空間上での目的関数に拡張される。 

\begin{align}
写像前の目的&関数 \\
& \underset{\alpha}{\max} \bigg[ -\dfrac{1}{2} \sum_{i=1}^n \sum_{j=1}^n \alpha_i \alpha_j y_i y_j \boldsymbol{x}_i^T \boldsymbol{x}_j + \sum_{i=1}^n \alpha_i \bigg] \\
& \\
特徴空間上に&拡張された目的関数 \\
& \underset{\alpha}{\max} \bigg[ -\dfrac{1}{2} \sum_{i=1}^n \sum_{j=1}^n \alpha_i \alpha_j y_i y_j \boldsymbol{\phi}(\boldsymbol{x}_i)^T \boldsymbol{\phi}(\boldsymbol{x}_j) + \sum_{i=1}^n \alpha_i \bigg] \\
\end{align}

この最適化問題を解けば、特徴空間上での分類境界を決定し、入力空間に戻すことで元データに対する非線形分離が可能となる。しかし、上の最適化問題をそのまま解くのは現実的には困難な場合がほとんどである。なぜなら、次元拡張により$~\boldsymbol{\phi}(\boldsymbol{x}_i)^T \boldsymbol{\phi}(\boldsymbol{x}_j)~$の内積部分の計算量が莫大になるためである。

カーネルトリック

次元拡張により計算量が膨大になった$~\boldsymbol{\phi}(\boldsymbol{x}_i)^T \boldsymbol{\phi}(\boldsymbol{x}_j)~$の内積部分の計算を簡略化するテクニックが__カーネルトリック__である。これは$~\boldsymbol{\phi}~$の内積部分を__カーネル関数__と呼ばれる関数で置き換える方法である。

K(\boldsymbol{x}_i,\boldsymbol{x}_j) = \boldsymbol{\phi}(\boldsymbol{x}_i)^T \boldsymbol{\phi}(\boldsymbol{x}_j)

このカーネル関数(内積なのでスカラー)を用いることで、2つの$~\boldsymbol{\phi}(\boldsymbol{x})~$(ベクトル)を直接計算することなく内積を見積もることが可能になり、特徴空間が高次元の場合でも双対問題を解く計算コストを大幅に削減することができる。

\underset{\alpha}{\max} \bigg[ -\dfrac{1}{2} \sum_{i=1}^n \sum_{j=1}^n \alpha_i \alpha_j y_i y_j K(\boldsymbol{x}_i,\boldsymbol{x}_j) + \sum_{i=1}^n \alpha_i \bigg], \quad \sum_{i=1}^n \alpha_i y_i = 0, \quad 0\leq\alpha_i ~~ (i=1, \cdots, n)

また、決定関数$~f(\boldsymbol{x})~$もカーネル関数を用いて

f(\boldsymbol{x})
 = \boldsymbol{w}^T\boldsymbol{\phi}(\boldsymbol{x}) + b
 = \sum_{i=1}^n \alpha_i y_i K(\boldsymbol{x}_i,\boldsymbol{x})

と表現することができる。つまり、カーネル関数さえ決めれば、最適化や決定関数の計算には$~\boldsymbol{\phi}(\boldsymbol{x})~$が一切必要なくなる。

カーネル関数

カーネル関数の代表的な関数形は、以下のの3つ。

関数に含まれる$~c, d, \gamma~$といったパラメータはハイパーパラメータで事前に決定しておく必要がある。また、ガウスカーネルのことは RBF (radial basis function) カーネルとよばれることもある。

\begin{align}
& 多項式カーネル\quad\quad\quad
  &K(\boldsymbol{x}_i,\boldsymbol{x}_j)
 &= [ \boldsymbol{x}_i^T\boldsymbol{x}_j + c ]^d \\
& ガウスカーネル
  &K(\boldsymbol{x}_i,\boldsymbol{x}_j)
 &= \exp ( -\gamma || \boldsymbol{x}_i - \boldsymbol{x}_j ||^2 ) \\
& シグモイドカーネル
  &K(\boldsymbol{x}_i,\boldsymbol{x}_j)
 &= \tanh ( b \boldsymbol{x}_i^T\boldsymbol{x}_j + c )
\end{align}

データに合わせて適切なカーネルを選んで計算を行うことで、SVM を用いた入力データの非線形分離が可能となります。

<参考動画>

実装演習

  • 各種モジュールをインポートする。
    スクリーンショット 2021-07-17 23.37.32.png

訓練データ生成① (線形分離可能)

  • 線形分離可能な訓練データを生成する。
    スクリーンショット 2021-07-17 23.37.53.png

  • 作成した訓練データを確認する。
    スクリーンショット 2021-07-17 23.50.46.png

学習

  • SVM で学習を行う。
    スクリーンショット 2021-07-17 23.39.14.png

予測

  • 変数を準備する。
    スクリーンショット 2021-07-17 23.39.33.png

  • 作られたモデルで予測を行う。
    スクリーンショット 2021-07-17 23.39.44.png

  • 予測結果をプロットする。
    スクリーンショット 2021-07-17 23.51.00.png

訓練データ生成② (線形分離不可能)

  • 線形分離可能な訓練データを生成する。
    スクリーンショット 2021-07-17 23.53.35.png

  • 作成した訓練データを確認する。
    スクリーンショット 2021-07-17 23.55.12.png

学習

  • SVM で学習を行う。
    スクリーンショット 2021-07-17 23.55.40.png

予測

  • 変数を準備する。
    スクリーンショット 2021-07-17 23.56.15.png

  • 作られたモデルで予測を行う。
    スクリーンショット 2021-07-17 23.56.27.png

  • 予測結果をプロットする。
    スクリーンショット 2021-07-17 23.56.39.png

ソフトマージンSVM

訓練データ生成③(重なりあり)

  • 重なりありの訓練データを生成する。
    スクリーンショット 2021-07-17 23.57.29.png

  • 作成した訓練データを確認する。
    スクリーンショット 2021-07-17 23.57.37.png

学習

  • SVM で学習を行う。
    スクリーンショット 2021-07-17 23.58.26.png

予測

  • 変数を準備する。
    スクリーンショット 2021-07-18 0.08.37.png

  • 作られたモデルで予測を行う。
    スクリーンショット 2021-07-18 0.08.45.png

  • 予測結果をプロットする。
    スクリーンショット 2021-07-18 0.08.53.png

0
1
0

Register as a new user and use Qiita more conveniently

  1. You get articles that match your needs
  2. You can efficiently read back useful information
  3. You can use dark theme
What you can do with signing up
0
1

Delete article

Deleted articles cannot be recovered.

Draft of this article would be also deleted.

Are you sure you want to delete this article?