線形回帰直線の求め方２~重回帰編~

はじめに

前回の記事では入力される変数と出力する変数がそれぞれ１つずつの単回帰を説明しました。今回は入力がn個あるときに求める重回帰直線を求めてみようと思ういます。とは言っても話はそこまで複雑になることはないと思います（行列の掛け算がわかれば難しくない…はず）。それではやっていきましょう。

重回帰直線の式を求める

前回の単回帰直線は、次のように書いていました

$$y = ax + b\tag{1}$$

ここではxが入力、yが出力で、求めたいパラメータはaとbでした。$$a=w_1, , b=w_0$$と文字の置き方を変えると、次のように書けますね。

$$y = w_0 + w_1x_1$$

ここで$$x_0=1$$として式を書き換えると、

$$y = w_0x_0 + w_1x_1\tag{2}$$

となりますね。ここまで来るとΣを使って表せそうです。こんなふうに。

$$y = \sum_{i=0}^{1}w_ix_i\tag{3}$$

これが入力が1のときの線形回帰の式です。では入力が2のときはどうなるかというと…

$$ y = w_0x_0 + w_1x_1 + w_2x_2 = \sum_{i=0}^{2}w_ix_i $$

じゃあ入力がnのときは…

$$y = \sum_{i=0}^{n}w_ix_i\tag{4}$$

となるわけです。Σの上部の数が変わっているだけですね。これで一般の線形回帰をあらわすことができました。あとは係数$$\bf w_i$$の求め方がわかれば予測を行うことができそうです。
先程の(4)式は行列を用いれば次のように表すことができます。

$$\bf{\hat{y}(w)} = \bf{\bar{x}} \bf{w}\tag{5}$$

また、単回帰のように係数は2乗誤差関数$$E(\bf{w})$$が最小化するようにするので、$$\frac{\partial E}{\partial\bf{w} }=0$$となるwを求めます。

$$E(\bf{w}) = || \bf{y} - \bf{\hat{y}(w)}||^2$$
$$=|| \bf{y} - \bf{\bar{x}w}||^2$$
$$= (\bf{y} - \bf{\bar{x}w})^T(\bf{y} - \bf{\bar{x}w})$$

ここで行列のノルムの2乗が行列の転置と行列の積で表されるのは、行列の形状を考えれば自明です。これによって$$E(\bf{w})$$はスカラになります。行列は形状を意識して計算をしなくてはならないので、例えば$$(\bf{\bar{x}w})^T =\bf{w^T\bar{x}^T} $$となることに注意します。式変形を続けます。

$$E(\bf{w}) = y^Ty - \bf{w^T\bar{x}^T}y - y^T\bar{x}w + \bf{w^T\bar{x}^T}\bar{x}w$$
$$= \bf y^Ty - 2y^Txw + ||\bar{x}w||^2$$
$$\nabla E = - 2\bf{y^T\bar{x}} + 2\bf {\bar{x}^T}\bar{x}w$$
$$(\bf{y^T\bar{x}})^T = \bar{x}^Ty$$
より、
$$\nabla E(\bf{w} )= -2\bf{\bar{x}^Ty} + 2\bar{x}^T\bar{x}w$$
$${E(\bf{w}}) = 0$$
とすれば良いので、
$$\bf{\bar{x}^T\bar{x}w }= \bf{\bar{x}^Ty} $$
両辺に$$\bf{(\bar{x}^T\bar{x})}^{-1}$$をかければ、次のように係数$\bf{w}$が求められます。
$${\bf{w} = \bf{(\bar{x}^T\bar{x})}^{-1}\bf{\bar{x}^Ty}}\tag{6}$$
これを実装していきます。

Pythonによる実装

#coding: utf-8  
import numpy as np  

class LinearRegression:  
    """  
    this class requires all input data's data type ndarray!!!  
    """  
    def __init__(self):  
        #係数wの初期化  
        self.w = None  

    def fit(self, x_train, y_train):  
        """  
        this function makes model fitting data.  
        """  
        x0 = np.ones((x_train.shape[0], 1))  
        x_train = np.concatenate([x0, x_train], axis=1)  
        self.w = np.linalg.inv(x_train.T.dot(x_train)).dot(x_train.T).dot(y_train)  

    def predict(self, x_test):  
        """  
        if you input data, this function predicts output values.  
        """  
        x0 = np.ones((x_test.shape[0], 1))  
        x_test = np.concatenate([x0, x_test], axis=1)  
        return x_test.dot(self.w)  

    def R2_score(self, x_test, y_test):  
        """  
        this function returns R2_score  
        """  
        x0 = np.ones((x_test.shape[0], 1))  
        x_test = np.concatenate([x0, x_test], axis=1)  
        pred = x_test.dot(self.w)  
        R2 = 1 - ( np.sum((y_test - pred)**2))/(np.sum((y_test - np.mean(y_test))**2))  
        return R2  

    def RMSE_score(self, x_test, y_test):  
        """  
        this function returns RMSE(Root Mean Squared Error) score  
        """  
        x0 = np.ones((x_test.shape[0], 1))  
        x_test = np.concatenate([x0, x_test], axis=1)  
        pred = x_test.dot(self.w)  
        RMSE = np.sum(np.sqrt((y_test - pred)**2)) / x_test.shape[0]  
        return RMSE

はじめて自作の関数作ったのであまりきれいなコードではないですが、とりあえず学習を行う関数fitと予測を行うpredictが大事です。R2_scoreはモデルの当てはまりの指標の決定係数を返し、RMSE_scoreは平均二乗誤差の平方根を返します。どちらも回帰モデルの評価によく使われます。
では次にこれがしっかり動くかをscikitlearnにあるデータセットを用いて確かめてみます。

Boston house price

米国ボストン市郊外における地域別の住宅価格のデータセットらしいです。各パラメータの解説はこちらが詳しいです。データセットの中身を確認するのはそちらを見たほうが良いと思います。今回は本筋から外れるので割愛して、とりあえず全部入力データとして突っ込んでみます。

from sklearn.datasets import load_boston  
from sklearn.model_selection import train_test_split  
boston = load_boston()  

x_train, x_test, y_train, y_test = train_test_split(boston.data, boston.target, random_state=0)  
model = LinearRegression.LinearRegression()  
model.fit(x_train, y_train)  
print(model.RMSE_score(x_test, y_test))    #3.669217386602595  
print(model.R2_score(x_test, y_test))      #0.6353620786674481

以上の結果から、RMSEは3.67、決定係数は0.635となりました。決定係数は1に近いほどデータと回帰の当てはまりが良いことになるので、このデータセットではただデータを突っ込むだけでは良いモデルを作ることは難しいようです。また、random_stateの値を変えればスコアは上下するので、これが最も良い線形回帰モデルというわけでもないです。訓練に使うデータによって大きく左右されています。