ラビットチャレンジ - 機械学習 第2章 非線形回帰モデル

# 0．概要
本記事は日本ディープラーニング協会認定の講座プログラムである「ラビット・チャレンジ」が提供している科目の1つである機械学習のレポートである。
記事タイトルに記載のとおり、第2章の非線形回帰モデルについて以下にまとめる。

# 1．非線形回帰モデル
データを線形(直線)で捉えられる場合は限られる。
複雑な非線形構造を持つ現象に対して行うモデリング。

## 1.1．基底展開法
回帰関数として、基底関数と呼ばれる既知の非線形関数とパラメータベクトルの線形結合を使用したもの。

```math
y_{i}=f(x_{i})+\epsilon_{i}
```

```math
y_{i}=w_{0}+\sum^{m}_{i=1}{w_{j}\phi_{j}(x_{i})}+\epsilon_{i}
```

$\phi_{j}(x_{i})$：基底関数
他の値は線形回帰と同じなので省略

以下が推定の流れ
1.説明変数
説明変数の数:$m$

```math
x_{i}=(x_{i1},x_{i2},...x_{i1m})\in\mathbb{R}^m
```

2.非線形関数ベクトル
基底関数の写像で変換し、$k$次元の特徴ベクトルとする。
基底関数の数:$k$

```math
\phi(x_{i})=(\phi_{1}(x_{i}),\phi_{2}(x_{i}),...\phi_{k}(x_{i}))^T\in\mathbb{R}^k
```

3.非線形回帰の計画行列

```math
\Phi^{(train)}=(\phi(x_{1}),\phi(x_{2}),...\phi(x_{n}))^T\in\mathbb{R}^{n{\times}k}
```

4.最尤法による予測
ここでは最尤法の例だが、最小二乗法でも同じように$\hat{y}$が求められる。

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

### 1.1.1．基底関数の例
* 多項式関数

```math
\phi_{j}=x^{j}
```

* ガウス型基底関数

```math
\phi_{j}(x)=exp\left(\frac{(x-\mu_{j})^2}{2h_{j}}\right)
```

2次元ガウス型基底関数の場合

```math
\phi_{j}(x)=exp\left(\frac{(x-\mu_{j})^{T}(x-\mu_{j})}{2h_{j}}\right)
```

* スプライン関数 / Bスプライン関数

# 2．未学習・過学習

## 2.1．未学習(アンダーフィッティング)
学習データに対して十分小さな誤差が得られないモデル(うまく学習できていない)。
例）線形回帰から非線形回帰に変更(データをプロットしてみたら曲線っぽかった等)

対策

* 表現力の高いモデルを利用する。

## 2.2．過学習(オーバーフィッティング)
学習データに対する誤差は小さくなったが検証データに対する誤差が大きいモデル(学習データの影響を受けすぎている)。
学習データに対しては精度良く予測できるが未知のデータに対して正しく予測できない。

対策

* 学習データを増やす。
* 不要な基底関数(説明変数)を削除して表現力を抑える。
* 正則化を使用して表現力を抑える。

# 3．正則化(罰則化法)
モデルの複雑さに伴って値が大きくなる正則化項を課した関数を最小化する。
ただ$MSE$等が最小となる点を探そうとすると過学習が起き得るため、原点近辺に制約を与えた上で最小となる点を探す。

正則化項がなければ最小二乗法等の推定値そのまま。

```math
S_{\gamma}=(y-\Phi{w})^T(y-\Phi{w})+\gamma{R}(w)
```

$S_{\gamma}$：正則化項を課した関数
$y$：目的変数
$\Phi$：計画行列
$w$：回帰係数
$\gamma$：正則化パラメータ。調整して決める。
$\gamma{R}(w)$：正則化項

## 3.1．L1ノルム
L1ノルムを利用した時の推定量をラッソ(Lasso)推定量と呼ぶ。
スパース推定とも呼ぶ。

いくつかのパラメータを0に推定する。
影響の小さい変数を0に推定するため、説明変数を選ぶ指標となる。


## 3.2．L2ノルム
L2ノルムを利用した時の推定量をリッジ(Ridge)推定量と呼ぶ。
縮小推定とも呼ぶ。

パラメータを0に近づけるよう推定する。
L1ノルムよりコストがかかるが精度向上に役立つ。

# 4．汎化性能
学習に使用したデータだけではなく、未知のデータに対する予測性能。
学習データとは別に用意した検証データで性能を測り、この検証データに対する誤差(汎化誤差)が小さいものが良い性能と言える。

## 4.1．訓練誤差
どれくらい正しく学習できたか(目的変数との差が小さいか)という指標。
小さくならないと学習できていないが、小さぎると過学習が起きている可能性あり。

```math
MSE_{train}=\frac{1}{n_{train}}\sum_{i=1}^{n_{train}}{(\hat{y}^{train}-y^{train})^2}
```

## 4.2．汎化誤差
未知のデータをどれだけ精度良く予測できたか(目的変数との差が小さいか)であり、モデルの性能を表す指標。

```math
MSE_{test}=\frac{1}{n_{test}}\sum_{i=1}^{n_{test}}{(\hat{y}^{test}-y^{test})^2}
```

汎化誤差は誤差、バイアス、バリアンスに分解できる。
低バリアンス・高バイアス：表現力が低いモデル
高バリアンス・低バイアス：表現力が高いモデル

### 4.2.1．バイアス
モデルの予測と目的変数との平均のズレ。

### 4.2.2．バリアンス
モデルの複雑さ(モデルの分散)。

## 4.3．ホールドアウト法
データを学習用とテスト用の2つ(学習用、検証用、テスト用に分けるなら3つ)に分割して使用する。
単純に分割するため、大量のデータがないと学習データが少なくてうまく学習できなかったり、テストデータが少なくて性能評価の精度が落ちたりする。
また、分割する時にデータの偏りがでないようにする。

線形回帰のハンズオンで7:3に分けていたのもこのホールドアウト法である。

## 4.4．交差検証法(クロスバリデージョン)
データをいくつかのブロックに分割し、1つのブロックをテストデータ、他を学習データとする(イテレータ)。
分割した数だけイテレータを繰り返してモデルをそれぞれ作っていくが、この時にテストデータとするブロックを順に変えていく。つまり、各イテレータの学習データ、テストデータは別々のものが使用されることになる。
できたモデルに対してそれぞれテストデータを与え、得られた精度の平均(CV値)をそれぞれ算出する。
各モデルのうち、最もCV値が小さいものが一番性能が良いと言える。

例）5分割した場合(n1～n5)
作られるモデルの数：5
イテレータの回数：5

イテレータ1：
　学習データ：n2～n5　テストデータ：n1

イテレータ2：
　学習データ：n1, n3～n5　テストデータ：n2

イテレータ3：
　学習データ：n1, n2, n4, n5　テストデータ：n3

イテレータ4：
　学習データ：n1～n3, n5　テストデータ：n4

イテレータ5：
　学習データ：n1～n4　テストデータ：n5

大元のデータは同じだが、それぞれのモデルは別々の学習データで学習していることから擬似的なアンサンブル学習である。

# 5．ハンズオン

## 5.1．ライブラリのインポート

```Python:
import numpy as np
import matplotlib.pyplot as plt
import seaborn as sns

%matplotlib inline
```

seabornはmatplotlibをよりきれいに表示したり簡易に実装できたりするラッパー的なライブラリ。

## 5.2．真の関数からノイズを伴うデータを生成

```Python:
n=100

def true_func(x):
    z = 1 - 48*x + 218*x**2 - 315*x**3 + 145*x**4
    return z 

def linear_func(x):
    z = x
    return z 

# 真の関数からノイズを伴うデータを生成

# 真の関数からデータ生成
data = np.random.rand(n).astype(np.float32)
data = np.sort(data)
target = true_func(data)

# ノイズを加える
noise = 0.5 * np.random.randn(n)
target = target + noise

# ノイズ付きデータを描画
plt.scatter(data, target)

plt.title('NonLinear Regression')
plt.legend(loc=2)
```

np.random.rand()を使ってn個生成した一様分布の乱数に対して真の関数を適用
真の関数：$1-48x+218x^2-315x^3+145x^4$

np.random.randn()を使ってn個生成した標準正規分布の乱数に対して0.5を掛けてノイズを生成

真の関数とノイズを加算した値ををプロットすると以下のようになる。
もちろん乱数を使っているため実行するたび結果が変わるが基本的な形は変わらない。

![nonlinear.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/389699/ee2c7c72-4334-c100-71a6-6d3a3a0f6e86.png)

## 5.3．線形回帰モデルで学習

```Python:
from sklearn.linear_model import LinearRegression

clf = LinearRegression()
data = data.reshape(-1, 1)
target = target.reshape(-1, 1)
clf.fit(data, target)

p_lin = clf.predict(data)

plt.scatter(data, target, label='data')
plt.plot(data, p_lin, color='darkorange', marker='', linestyle='-', linewidth=1, markersize=6, label='linear regression')
plt.legend()
print(clf.score(data, target))
```

> 0.4501186554569356

![nonlinear2.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/389699/efd885b4-932f-fc64-b2f7-bbc267b37847.png)

データを見ると明らかに曲線的であることからも分かる通り、線形回帰モデルは適していない。また、決定係数も0.45程度となっている。

reshapeはnumpy配列の変形をする関数だが、-1を指定すると元の形状から推測された形状に変換される。

```Python:
data2 = np.random.rand(n).astype(np.float32)
data2 = np.sort(data2)
target2 = true_func(data2)

print(data2.shape)
print(target2.shape)

print(data2.reshape(-1, 1).shape)
print(target2.reshape(-1, 1).shape)
```

> (100,)
> (100,)
> (100, 1)
> (100, 1)

変数の値が変わると後続のグラフも変わってしまうので別変数で形状を見てみるとreshapeする前は1次元だと分かる。
reshape(-1, 1)をすることで推測されたサイズ(データ数)×1に変形される。

結果的になにも変わっていないようにも思えるがreshapeすることで2次元配列となっている。
fit関数へ与えるデータは1次元ではエラーとなるためこれで実行できる形に変形している。

## 5.4．カーネルリッジ回帰モデルで学習(KernelRidgeを使用)

```Python:
from sklearn.kernel_ridge import KernelRidge

clf = KernelRidge(alpha=0.0002, kernel='rbf')
clf.fit(data, target)

p_kridge = clf.predict(data)

plt.scatter(data, target, color='blue', label='data')

plt.plot(data, p_kridge, color='orange', linestyle='-', linewidth=3, markersize=6, label='kernel ridge')
plt.legend()

print(clf.score(data, target))
```

> 0.8625198512970995

![nonlinear3.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/389699/5f2b5f5e-659a-b734-8962-6fc841f41159.png)

カーネルリッジ回帰モデルとはカーネル回帰にL2ノルム(リッジ推定量)を追加したモデル。
カーネル回帰は非線形回帰モデルの1つ。
グラフを見ると分かる通り、良い感じに適合している。決定係数も0.86と線形回帰の2倍近くもある。

参考
https://qiita.com/wsuzume/items/09a59036c8944fd563ff
https://masamunetogetoge.com/kernelpart1
https://masamunetogetoge.com/overfit-and-generalization

alpha：正則化パラメータ
kernel：カーネル関数の種類。

## 5.5．カーネルリッジ回帰モデルで学習(KernelRidgeを不使用)

```Python:
from sklearn.metrics.pairwise import rbf_kernel
from sklearn.linear_model import Ridge

kx = rbf_kernel(X=data, Y=data, gamma=50)

clf = Ridge(alpha=30)
clf.fit(kx, target)

p_ridge = clf.predict(kx)

plt.scatter(data, target,label='data')
for i in range(len(kx)):
    plt.plot(data, kx[i], color='black', linestyle='-', linewidth=1, markersize=3, label='rbf', alpha=0.2)

plt.plot(data, p_ridge, color='green', linestyle='-', linewidth=1, markersize=3,label='ridge regression')

print(clf.score(kx, target))
```

> 0.8483107400405749

![nonlinear4.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/389699/af3f29c3-a55a-3916-294c-f261a814233e.png)

rbf_kernelに説明変数を与え、RBFカーネルと目的変数をリッジ回帰に与えている。
ハイパーパラメータの違いによるものだと思われるがKernelRidgeを使った時より多少決定係数の値は劣る。

gamma：カーネルの形らしい

参考：https://qiita.com/kotaroito/items/4eb29d42d7f8c534332f

## 5.6．多項式回帰モデルで学習

```Python:
deg = [1, 2, 3, 4, 5, 6, 7, 8, 9, 10]

for d in deg:
    regr = Pipeline([
        ('poly', PolynomialFeatures(degree=d)),
        ('linear', LinearRegression())
    ])
    regr.fit(data, target)
    
    # make predictions
    p_poly = regr.predict(data)
    
    # plot regression result
    plt.scatter(data, target, label='data')
    plt.plot(data, p_poly, label='polynomial of degree %d' % (d))
    print(regr.score(data, target))
```

> 0.450118667201583
> 0.5836621504744504
> 0.6398637941304797
> 0.8651962282474446
> 0.8683762396358908
> 0.8697061688538484
> 0.8756350078101446
> 0.8757686555122806
> 0.8765569664601603
> 0.875972484556606

![nonlinear5.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/389699/b6905bdd-06f4-6c45-eeb1-2f2ea532e92e.png)


PolynomialFeaturesは指定した次数を元に説明変数を調整してくれる。
モデルクラスとしてはLinearRegressionを使っているので正則化項はなし。

1～10次それぞれでじっこうしているため決定係数もそれぞれ出ている。
5.3でやったとおり、1次では全然適合できていないが3次以降はあまり変わっていないことが分かる。
むしろ10次は9次より若干劣っている。

## 5.7．カーネル回帰+ラッソ回帰

```Python:
from sklearn.metrics.pairwise import rbf_kernel
from sklearn.linear_model import Lasso

kx = rbf_kernel(X=data, Y=data, gamma=5)

lasso_clf = Lasso(alpha=10000, max_iter=1000)
lasso_clf.fit(kx, target)

p_lasso = lasso_clf.predict(kx)

plt.scatter(data, target)

plt.plot(data, p_lasso, color='green', linestyle='-', linewidth=3, markersize=3)

print(lasso_clf.score(kx, target))
```

> -2.220446049250313e-16

![nonlinear6.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/389699/00a5ea9e-238f-0c0e-387e-b7fb91d0a4ab.png)

KernelRidgeを使わずにカーネルリッジ回帰を行った5.5のラッソ回帰版。
グラフを見れば一目瞭然だが決定係数もマイナスになるほど全く適合していない。

ググってもカーネルラッソ回帰という言葉はひっかからなかったのでそんなものはない？

## 5.8．サポートベクター回帰モデル(SVR)

```Python:
from sklearn import model_selection, preprocessing, linear_model, svm

# SVR-rbf
clf_svr = svm.SVR(kernel='rbf', C=1e3, gamma=0.1, epsilon=0.1)
clf_svr.fit(data, target)
y_rbf = clf_svr.fit(data, target).predict(data)

print(clf_svr.score(data, target))
 
# plot
plt.scatter(data, target, color='darkorange', label='data')
plt.plot(data, y_rbf, color='red', label='Support Vector Regression (RBF)')
plt.legend()
plt.show()
```

![nonlinear7.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/389699/b77e0604-3af4-cbf3-57a2-06b24a18ed56.png)

サポートベクターマシンの回帰版。
おおよその形には沿っているがあまり適合していない。

C：ペナルティ項(正則化項？)の係数
epsilon：どれだけ回帰直線と離れていたらペナルティをかけるか

参考
http://neuro-educator.com/ml19/

## 5.9．Kerasによる深層学習実装

```Python:
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/study_ai_ml/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/study_ai_ml/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))

# 予測
y_pred = estimator.predict(x_train)
print(y_pred)

plt.title('NonLiner Regressions via DL by ReLU')
plt.plot(data, target, 'o')
plt.plot(data, true_func(data), '.')
plt.plot(x_train, y_pred, "o", label='predicted: deep learning')

print(lasso_clf.coef_)
```

y_pred
> 18/18 [==============================] - 0s 1ms/step
> [ 1.2567639   0.36203647  0.45334932 -2.123969   -2.015928   -1.4220177
>  -2.2667038  -2.2294168  -1.3345145  -2.025425   -0.34206712  1.0893227
>  -0.86253875  0.84631246  0.9902647   1.210631   -1.8815542   0.38569453
>   0.35375044  0.51373094  1.1206777   0.37869334  0.35905007  0.93582046
>   1.2273624  -2.1871383  -2.2800443   1.2619673   1.0792997   0.6259613
>  -0.5350753   1.0917804  -0.3922692  -0.422297   -1.1384251  -2.170742
>  -2.2828817  -0.49975392 -0.5413679  -1.3219417   0.35937083  1.1627012
>  -0.5350753  -1.0982922  -2.1204653   1.2239919  -2.0876467  -1.0516531
>  -2.2573252   1.234646   -2.0866053  -1.4121388   0.9806428   1.242071
>   0.4051356   0.35369515 -1.3731359   0.35273564 -2.1285925   0.35879287
>  -2.2795181  -1.9124163   1.2147804   0.38208917  0.3823473   0.9335167
>  -2.2768147   0.36280116  0.8141983   0.46964574 -2.228463   -0.5908961
>   0.37028825  1.1164801  -0.51945305  1.0124962   1.2267804  -2.123545
>   0.35783193  0.8735486   0.35455173 -0.7707679  -2.268309    0.35305035
>  -2.0749183   0.38055658  0.96130353  1.0345464   1.1318389   0.9887817 ]

coef_
> [-0. -0. -0. -0. -0. -0. -0. -0. -0. -0. -0. -0. -0. -0. -0. -0. -0. -0.
>  -0. -0. -0. -0. -0. -0. -0. -0. -0. -0. -0. -0. -0. -0. -0. -0. -0. -0.
>  -0. -0. -0. -0. -0. -0. -0. -0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.
>   0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.
>   0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.
>   0.  0.  0.  0.  0.  0.  0.  0.  0.  0.]

![nonlinear8.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/389699/f2354526-d3be-a86b-9637-25cbb6fe1e4a.png)

KerasはTensorflow等のラッパーライブラリ。

レイヤーを追加していくようにニューラルネットワークを組んだりできるライブラリ。
入力層は全結合層ユニット数10で活性化関数が正規化線形関数(ランプ関数)。
隠れ層として全結合層ユニット数1000で活性化関数が正規化線形関数(ランプ関数)が7層+全結合層ユニット数1000で活性化関数が線形関数の層が1つ。
出力層としてユニット数1つの全結合層。

損失関数としてMSE。最適化関数としてAdam。

# X．ラビットチャレンジとは
ラビットチャレンジとは、日本ディープラーニング協会認定の講座プログラムの1つ。
E資格を受験するためにはこのラビットチャレンジ等、いずれかの講座プログラムを修了しなければならない。

ラビットチャレンジの特徴は「現場で潰しが効くディープラーニング講座」の通学講座録画ビデオを編集した教材を使用した自習スタイルであるという点。
サポートは他の講座より少なく、受け身ではなく自主的に学んでいく姿勢でなければ進められないが、その分、他の講座に比べると安価であり、手が出しやすい。
ある程度知識がある人、自力で頑張るぞというガッツのある人向けではないかと感じる。

![bnr_jdla.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/389699/bce66671-a2f5-a7f9-de6f-3f87e15c253b.png)