回帰
（線形）単回帰モデル
1種類のデータから1種類の予測値を出力するモデルです。
入力データを$$x=(x_1, x_2, \cdots, x_n)$$、出力データを$$y=(y_1, y_2, \cdots, y_n)$$とするとき、モデルを最もよく表している__直線$$y=wx+b$$の傾き$w$と切片$b$を決めます。
$w$を__重み、$b$を__バイアス__といいます。
実際のデータとモデル（$y=wx+b$）の間には誤差があるので、その誤差を
$$\varepsilon=(\varepsilon_1, \varepsilon_2,\cdots , \varepsilon_n)$$とします。このとき、各教師データについて
y_i = wx_i+b+\varepsilon_i
が成り立ちます。
モデルを最もよく表しているとは、適切な__損失関数__（例えば、二乗和誤差）を決め、モデル（直線の式）がそれを最小にするようなときにいいます。
損失関数$L(w,b)$を二乗和誤差（残差平方和ともいう）
L(w, b):=\sum_{i=1}^n \{y_i-(wx_i+b)\}^2
にすると、$$L(w,b)=\sum_{i=1}^n\varepsilon_i^2$$すなわち誤差の二乗和となっていることが分かります。場合によっては、(1/2)倍します（微分すると係数が1にできるため）。
詳細は省きますが、$L(w,b)$を$w, b$それぞれで偏微分した式$=0$を解くことで、$w, b$を求めることができます。
b = \bar{y}-w\bar{x}, \ w=\frac{Cov(x, y)}{Var(x)}=\frac{\sum_{i=1}^n(x_i-\bar{x})(y_i-\bar{y})}{\sum_{i=1}^n(x_i-\bar{x})}
ここで、$\bar{x}$は$x$の平均、$Var(x)$は$x$の分散、$Cov(x, y)$は$x$と$y$の共分散を表しています。
特に、この$(w, b)$の元で
\begin{eqnarray}
&&y=wx+b=wx+(\bar{y}-w\bar{x})=\bar{y}+w(x-\bar{x})\\
&\Leftrightarrow& y-\bar{y}=w(x-\bar{x})
\end{eqnarray}
より、回帰直線は観測データの標本平均$(\bar{x}, \bar{y})$を通ります。
実装
Pythonで自力で実装＆sklearnのlinear_modelのLinearRegressionモデルを用いた実装をしてみます。
なお、再現性のために乱数シードを設定しています。
まずは必要なライブラリをインポート。
import numpy as np
import matplotlib.pyplot as plt
from numpy.random import randn
次に、線形回帰を関数として定義。
def LinearRegression(x, y):
n = len(x)
temp_x = x.sum()/n
temp_y = y.sum()/n
w = ((x-temp_x)*(y-temp_y)).sum()/((x-temp_x)**2).sum()
b = temp_y - w * temp_x
return w, b
そして、教師データとして、$y=x+1+noise$を想定。$noise$は適当に標準正規分布に従う乱数に0.1をかけたものとしました。
# initial value
# y=x+1にノイズを加えたものを想定
x = np.array([1,2,3,4,5])
np.random.seed(seed=0)
y = np.array([2,3,4,5,6])+np.random.randn(5)*0.1
print(y)
w, b = LinearRegression(x, y)
print(w, b)
これを実行すると、およそ$w=1.02, b=1.08$と求まります。
念のため、グラフで可視化します。
plt.scatter(x, y, color="k")
xx = np.linspace(0, 5, 50)
yy = xx * w + b
plt.plot(xx, yy)
plt.show()
いい感じですね！
ちょっとコードを使いやすいように書き換える
クラスを使ってコードを書き換えて見ます。
modelというクラスを作り、パラメータを求めるfitメソッド、予測を出力するpredictメソッド、決定係数を計算するscoreメソッドを作りました。決定係数$R^2$は、1に近いほど精度がよいということです。
R^2:= 1-\frac{\sum_{i=1}^n (y_i-(wx_i+b))^2}{\sum_{i=1}^n (y_i-\bar{y})^2}
と定義されます。
class model():
def fit(self, X, y):
n = len(X)
temp_X = X.sum()/n
temp_y = y.sum()/n
self.w = ((X-temp_X)*(y-temp_y)).sum()/((X-temp_X)**2).sum()
self.b = temp_y - self.w * temp_X
def predict(self, X):
return self.w*X +self.b
def score(self, X, y):
n = len(X)
return 1 - ((y-(self.w*X+self.b))**2).sum() / ((y-y.sum()/n)**2).sum()
先ほどと同じように、教師データを入力します。
加えて、テストもさせてみます。
X = np.array([1,2,3,4,5])
np.random.seed(seed=0)
y = np.array([2,3,4,5,6])+np.random.randn(5)*0.1
lr = model()
lr.fit(X, y)
print("w, b={}, {}".format(w, b))
# テストデータ
test_X = np.array([6,7,8,9,10])
pre = lr.predict(test_X)
print(pre)
# 決定係数R^2
print(lr.score(X, y))
グラフ化します。
# 結果をグラフに表示
plt.scatter(X, y, color="k")
plt.scatter(test_X, pre, color="r")
xx = np.linspace(0, 10, 50)
yy = xx * w + b
plt.plot(xx, yy)
plt.show()
sklearnのライブラリを利用する
最後に、linear_modelの線形回帰LinearRegression()を用いて簡単にコードを書いてみましょう。
from sklearn import linear_model
X = np.array([1,2,3,4,5])
np.random.seed(seed=0)
y = np.array([2,3,4,5,6])+np.random.randn(5)*0.1
X = X.reshape(-1,1)
model = linear_model.LinearRegression()
model.fit(X, y)
print(model.coef_[0])
print(model.intercept_)
# 決定係数R^2
print(model.score(X, y))
$w$はcoef_[0]、$b$はintercept_に相当します。値が一致していることがわかると思います。ここで注意する点は、1次元の入力データの場合は$X=np.array([[1],[2],[3],[4],[5]])$のように入力する必要があります。$X=np.array([1,2,3,4,5])$だとエラーを吐きます。
次に、先ほどのテストデータで予測をしてみます。
test_X = test_X.reshape(-1,1)
pre = model.predict(test_X)
グラフ化します。
plt.scatter(X, y)
plt.scatter(test_X, pre)
plt.plot(xx, model.coef_[0]*xx+model.intercept_)
plt.show()
まとめ
sklearnを用いてコードを書くと、簡単に次のようにかけます。これは、他の手法でもほぼ同様の形でかけます。あとはパラメータであったり、教師データを訓練データとテストデータに分けるなどオプション的な要素が加わるだけです。
# モデルを定義
model = linear_model.LinearRegression()
# 学習(X:入力データ、y:正解ラベル)
model.fit(X, y)
# 予測(XX:未知のデータ)
model.predict(XX)
単回帰（最小二乗法）は以上です。
次回に続く。。。
追記
よく考えると、$b=(b, b, \cdots, b)\in \mathbf{\mathrm{R}}^n$とすべきでしたね。Pythonコードではブロードキャストのおかげでスカラー$b$としてもベクトルや行列との和は自動的に計算できます。
追記その2
from sklearn import linear_model
としてますが、
from sklearn.linear_model import LinearRegression
の方が楽ですね。他の書籍でもこちらが一般的かと。
参考文献
- 加藤公一『機械学習のエッセンス』SB Creative、2019年