scikit-learnで線形モデルとカーネルモデルの回帰分析をやってみた - イラストで学ぶ機会学習

[Amazon: イラストで学ぶ 機械学習 最小二乗法による識別モデル学習を中心に](https://www.amazon.co.jp/%E3%82%A4%E3%83%A9%E3%82%B9%E3%83%88%E3%81%A7%E5%AD%A6%E3%81%B6-%E6%A9%9F%E6%A2%B0%E5%AD%A6%E7%BF%92-%E6%9C%80%E5%B0%8F%E4%BA%8C%E4%B9%97%E6%B3%95%E3%81%AB%E3%82%88%E3%82%8B%E8%AD%98%E5%88%A5%E3%83%A2%E3%83%87%E3%83%AB%E5%AD%A6%E7%BF%92%E3%82%92%E4%B8%AD%E5%BF%83%E3%81%AB-KS%E6%83%85%E5%A0%B1%E7%A7%91%E5%AD%A6%E5%B0%82%E9%96%80%E6%9B%B8-%E6%9D%89%E5%B1%B1/dp/4061538217)を読み進めています。

![イラストで学ぶ機械学習](https://images-fe.ssl-images-amazon.com/images/I/6106R-NIuzL._AC_UL115_.jpg "イラストで学ぶ機械学習")

出典： [Amazon](https://images-na.ssl-images-amazon.com/images/I/6106R-NIuzL._SX355_BO1,204,203,200_.jpg)



今回は第３章〜第5章の最小二乗学習による回帰分析をpythonとsckit-learnで実装してみます。利用環境は以下の通りです。

- anaconda3-2.5.0
- scikit-learn (0.17.1)


## データを用意する

右肩上がりのトレンド（第二項）に周期変動（第一項）とノイズ（第三項）が乗っているというデータを用意します。本に比べてノイズパラメタを増やしています。

![sin.png](https://qiita-image-store.s3.amazonaws.com/0/20304/8895e846-d417-e067-ef7b-60c8359d7219.png)

```py3
import numpy as np
import matplotlib.pyplot as plt

%matplotlib inline

n = 50; N = 1000

x = np.linspace(-3, 3, n)
X = np.linspace(-3, 3, N)

pix = np.pi * x
y = np.sin(pix) / pix + 0.1 * x + 0.2 * np.random.randn(n)

x = x.reshape(-1, 1)
X = X.reshape(-1, 1)
y = y.reshape(-1, 1)

plt.scatter(x,y)
```

## LinearRegression

何も制約が付かないプレーンな線形モデルである[sklearn.linear_model.LinearRegression](http://scikit-learn.org/stable/modules/generated/sklearn.linear_model.LinearRegression.html)を利用して、フィッティングしてみます。

![linear.png](https://qiita-image-store.s3.amazonaws.com/0/20304/963afdfb-f9b1-547c-cd6d-00477936c1e9.png)


```py3
from sklearn.linear_model import LinearRegression

clf = LinearRegression()
clf.fit(x, y)

p = clf.predict(X)

plt.scatter(x, y)
plt.plot(X,p)

print(clf.score(x, y))
```

なんとなく右肩上がりのトレンドは見えますが、決定係数(coefficient of determination)は0.10と低いです。



## KernelRidge

周期的に変動するデータは線形モデルと相性が悪いので、次はガウスカーネルモデルを利用します。

お手軽にやるなら[sklearn.kernel_ridge.KernelRidge](http://scikit-learn.org/stable/modules/generated/sklearn.kernel_ridge.KernelRidge.html) がよさそうです。

> Kernel ridge regression (KRR) combines ridge regression (linear least squares with l2-norm regularization) with the kernel trick. It thus learns a linear function in the space induced by the respective kernel and the data. For non-linear kernels, this corresponds to a non-linear function in the original space.

英文の説明そのままですが、Kernel Ridge Regression は l2制約付き最小二乗学習（linear least squares with l2-norm regularization）です。

![kernel_ridge.png](https://qiita-image-store.s3.amazonaws.com/0/20304/abe2d2f0-ed87-2686-11e2-8569838c0bdc.png)


```py3
from sklearn.kernel_ridge import KernelRidge

clf = KernelRidge(alpha=1.0, kernel='rbf')
clf.fit(x, y)

p = clf.predict(X)

plt.scatter(x, y)
plt.plot(X, p)

print(clf.score(x, y))
```

決定係数は0.81で、線形モデルから大幅に改善しました。
余談ですが、ドキュメントがなかなか見つからずに苦労しました ^^

https://github.com/scikit-learn/scikit-learn/blob/51a765a/sklearn/kernel_ridge.py#L120-L121 を追っていくと、kernelに'rfb'を指定できることがわかります。


## KernelRidgeを使わない最小二乗学習

KernelRidgeを使わずに、RBFカーネルによる最小二乗学習を実行することもできます。


### sklearn.linear_model.LinearRegression

まず制約を何も付けずにやってみます。

```py3
from sklearn.metrics.pairwise import rbf_kernel

kx = rbf_kernel(x, x)
KX = rbf_kernel(X, x)

clf = LinearRegression()
clf.fit(kx, y)

p = clf.predict(KX)

plt.scatter(x, y)
plt.plot(X, p)

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

![linear_regression_kernel.png](https://qiita-image-store.s3.amazonaws.com/0/20304/b009c79c-766a-afd3-c908-cb0eec4e51c7.png)

オーバーフィット気味です。

### sklearn.linear_model.Ridge

次に、l2制約を入れてみます。

```py3
from sklearn.metrics.pairwise import rbf_kernel
from sklearn.linear_model import Ridge

kx = rbf_kernel(x, x)
KX = rbf_kernel(X, x)

clf = Ridge()
clf.fit(kx, y)

p = clf.predict(KX)

plt.scatter(x, y)
plt.plot(X, p)

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

![ridge_with_kernel.png](https://qiita-image-store.s3.amazonaws.com/0/20304/efd67676-6c67-85d4-a7f7-12b0062c602d.png)

滑らかになりました。


### sklearn.linear_model.Lasso

今度はl1制約付き（ラッソ回帰）でやってみます。

```py3
from sklearn.metrics.pairwise import rbf_kernel
from sklearn.linear_model import Lasso

kx = rbf_kernel(x, x)
KX = rbf_kernel(X, x)

clf = Lasso(alpha=0.01)
clf.fit(kx, y)

p = clf.predict(KX)

plt.scatter(x, y)
plt.plot(X, p)

print(clf.score(kx, y)) # 0.820550922167
```

![lasso_with_rbf_kernel.png](https://qiita-image-store.s3.amazonaws.com/0/20304/37a59066-dbfc-7228-ef86-88756bbebd66.png)


ラッソ回帰はパラメタをスパースにするので、確認してみると...

```py3
print(clf.coef_) 

[-0.         -0.         -0.         -0.         -0.         -0.         -0.
 -0.         -0.         -0.         -0.         -0.         -0.         -0.
 -0.05432222 -0.         -0.         -0.         -0.         -0.          0.
  0.          0.          0.          0.17642978  1.08522759  0.          0.
  0.          0.          0.         -0.         -0.         -0.         -0.
 -0.         -0.         -0.         -0.          0.          0.          0.
  0.          0.          0.          0.          0.66599521  0.01403876
  0.          0.        ]
```

期待通り 0 が多くなってます。


## sklearn.linear_model.ElasticNet

モデルを挿げ替えるだけですが、一応l1 + l2制約付き（ElasticNet）でやってみます。

```py3
from sklearn.linear_model import ElasticNet

kx = rbf_kernel(x, x)
KX = rbf_kernel(X, x)

clf = ElasticNet(alpha=0.01)
clf.fit(kx, y)

p = clf.predict(KX)

plt.scatter(x, y)
plt.plot(X, p)

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

![elastic_net.png](https://qiita-image-store.s3.amazonaws.com/0/20304/c3038c69-d286-5253-0722-80eb6bd9e7a5.png)

今回の例では目視で確認できるほどの差異はないようです。

## まとめ

周期変動するダミーデータを利用して、LinearRegression、Ridge、Lasso、ElasticNet など様々な回帰分析をやってみました。