最小2乗法(OLS)をscikit-learnで使ってみる

今回は、回帰モデルで最も基本的なもの線形回帰モデルをscikit-learnで使うことにする。
以下の公式ページを参考にしている。

1.1 Generalized Linear Models – scikit-learn

最小2乗回帰(Ordinary Least Squares)
- 定式化
- 実装
まとめ

最小2乗回帰(Ordinary Least Squares)

定式化

与えられた
$$
x_1, x_2, \cdots, x_p \
$$
というp個の値が与えられた状態で値yを予測するモデル
$$
\hat{y}(w,x) = w_o + w_1x_1 + \cdots + w_px_p \
$$ を構築することが目的。
なので、学習させたいのはパラメータ群
$$
w = (w_0, w_1, \cdots, w_p) \
$$
となる。
最小2乗回帰では、モデルの予測値と正解値yの差が最小になるようなwの値の組を求めればよい。
定式化すると求めるwを
$$
\hat{w} \
$$
とすると
$$
\hat{w} = \underset{w}{argmin} |Xw – y|^2
$$
ここでXとyは
$$
y = (y_1, y_2, \cdots y_N)^T \
x_i = (x_{i0}, x_{i1}, x_{i2}, \cdots, x_{ip})^T \
(x_{i0} = 1) \
X = (x_1^T, x_2^T, \cdots, x_N^T)^T \
$$
となり、
$$
x_iとy_iがそれぞれ対応している　　
$$
このとき、wについて偏微分した値がゼロになればよいので、
$$
\frac{\partial}{\partial w} |Xw – y | ^2 = 0 \
w = (X^T X)^{-1} X^T y \
$$
これは
$$
X^TX
$$
に逆行列が存在するときにのみ存在する。

実装

とりあえずこれをscikit-learnで使えるかどうか確かめてみる。
以下のようなモデルを学習できるかどうか確かめる。
$$
y = -3.0 + 2.0x_1 + noise \
$$
ノイズ項を入れることで予測を若干難しくする。
まずは学習する元のデータを作成。


import numpy as np
np.random.seed(42)
w = np.array([2.0])
x = np.random.randn(10).reshape((10,1))
noise = np.random.randn(10)/2.0 # 標準偏差を半分程度に抑える
noise = noise.reshape((10,1))
y = x*w -3.0 + noise

これをグラフにプロットする。


import matplotlib.pyplot as plt
plt.scatter(x, y)
plt.xlabel('x')
plt.ylabel('y')
plt.grid()
plt.show()

これにフィッティングするような直線を学習させる。
ここではscikit-learnのLinearRegressionを使用する。


from sklearn import linear_model
reg = linear_model.LinearRegression()  
reg_fit = reg.fit(x, y)

学習したパラメータを見てみる。


print("coefficient is \n", reg.coef_)
print("intercept is \n", reg.intercept_)

結果


coefficient is
 [[1.94227475]]
intercept is
 [-3.36946468]

ほどほどに学習していそう。(interceptがw_0、coefficientがw_1に該当)
グラフを重ねてプロットしてみる。


import matplotlib.pyplot as plt
x_lin = np.linspace(-0.5, 1.5, 20).reshape(-1,1)
plt.scatter(x, y, label='raw data')
plt.plot(x_lin, reg_fit.predict(x_lin), label='predict')
plt.xlabel('x')
plt.ylabel('y')
plt.legend()
plt.grid()
plt.savefig('ols_predict.png')
plt.show()

一変量だとこんな感じ。
多変量だと以下のようになる。
$$
y = -1.5 + 2.0x_1 – 3.0x_2 \
$$
を学習させるとする。


np.random.seed(42)
w_2 = np.array([2.0, -3.0])
x_2 = np.random.randn(15,2)  
noise = np.random.randn(15)/2.0  
y_2 = np.dot(x_2, w_2) -1.5 + noise

これを学習させる。


reg_2 = linear_model.LinearRegression().fit(x_2, y_2)
print("w_0 is \n", reg_2.intercept_)
print("w_1 w_2 is \n", reg_2.coef_)

結果。


w_0 is
 -1.587149889988845
w_1 w_2 is
 [ 2.22325845 -2.94948391]

数式
$$
w = (X^T X)^{-1} X^T y \
$$
で計算してみると、


X = np.ones((15,3))  
X[:, 1:] = x_2
w_pred = np.dot(np.dot(np.linalg.inv(np.dot(X.T, X)), X.T), y_2.reshape(-1,1))
print(w_pred)

結果


[[-1.58714989]
 [ 2.22325845]
 [-2.94948391]]

計算結果がほぼ同じになったことが分かる。

まとめ

今回は線形回帰モデルの中でも最も基礎的な最小二乗法について扱った。
scikit-learnだと数式がよく分かってなくてもとりあえず計算できるのがありがたい。(このあたりはエクセルでも簡単にできるものではあるが)
次は罰則項として設定したものに係数の大きさで制限をかけるタイプのモデル(LassoとRidge)について扱っていく。

Pocket