【Python】重回帰分析を最小二乗法で実装する方法を解説

Pythonを使って、重回帰分析を自分で実装してみたい…。

今回はこんな人のために、重回帰分析を最小二乗法で実装する方法を解説しています。

本記事の流れは次の通りです。

＊最小二乗法や線形代数の基礎は理解している、という前提で解説しているので予めご了承ください。

最小二乗法で係数を求める【線形代数】

\(x_1 ^{(i)} ,x_2 ^{(i)} ,\cdots ,x_k ^{(i)} \)：i番目のデータの説明変数の値
\(x_0^{(i)}\)：すべての\(i\)で1（計算の都合上）
\(w_1,w_2,\cdots w_k\)：回帰係数
\(w_0\)：切片

損失関数を行列で表す

最小二乗法は「正解の値」と「回帰モデルで予測した値」の差の平方和を最小化する、というものでした。

i番目のデータの正解の値を\(y_i\)，予測値を\(\hat{y_i}=w_0x_0^{(i)}+w_1x_1 ^{(i)} +w_2x_2 ^{(i)} +\cdots +w_kx_k ^{(i)}\)とすると，損失関数(Loss関数，コスト関数,評価関数などともいう)Lは

$$L=\sum_{i=1}^{n}(y_i-\hat{y_i})^2$$

となります。これを最小化してくれるような回帰係数 \(w_1,w_2,\cdots w_k\) を求めるために，行列を用いて\(L\)をキレイな形に書きなおしたいと思います。

まず，

$$\boldsymbol{y} = \begin{pmatrix}y_0\\y_1\\\vdots\\y_n\end{pmatrix} \ \ ,\ \ \boldsymbol{w} = \begin{pmatrix}w_0\\w_1\\\vdots\\w_k\end{pmatrix}$$
$$\boldsymbol{X} = \begin{pmatrix}x_0^{(1)} & x_1^{(1)} & x_2^{(1)} & \cdots & x_k^{(1)} \\ x_0^{(2)} & x_1^{(2)} & x_2^{(2)} & \cdots & x_k^{(2)}\\ \vdots&\vdots&\vdots&\ddots&\vdots\\x_0^{(n)} & x_1^{(n)} & x_2^{(n)} & \cdots & x_k^{(n)} \end{pmatrix} $$

とおくと，損失関数\(L\)は次のように表せます。

$$\begin{align*} L &= \begin{pmatrix}y_1-\hat{y_1}&y_2-\hat{y_2}&\cdots &y_n-\hat{y_n}\end{pmatrix} \begin{pmatrix} y_1-\hat{y_1}\\y_2-\hat{y_2}\\\vdots \\y_n-\hat{y_n} \end{pmatrix}\\ &= (\boldsymbol{y}- \boldsymbol{\hat{y}})^{T}(\boldsymbol{y}- \boldsymbol{\hat{y}}) \end{align*}$$

回帰モデルが\(\boldsymbol{\hat{y}}=\boldsymbol{Xw}\)となることから，

$$ \begin{align*}L &= (\boldsymbol{y}- \boldsymbol{\hat{y}})^{T}(\boldsymbol{y}- \boldsymbol{\hat{y}}) \\
&= (\boldsymbol{y}- \boldsymbol{Xw} )^{T}(\boldsymbol{y}- \boldsymbol{Xw}) \\
&= (\boldsymbol{y}^T – (\boldsymbol{Xw})^T )(\boldsymbol{y}- \boldsymbol{Xw}) \\ &= (\boldsymbol{y}^T – \boldsymbol{w}^T \boldsymbol{X}^T )(\boldsymbol{y}- \boldsymbol{Xw}) \\
&= \boldsymbol{y}^T \boldsymbol{y} \ -\ \boldsymbol{y}^T \boldsymbol{Xw}\ -\ \boldsymbol{w}^T\boldsymbol{X}^T\boldsymbol{y}\ +\ \boldsymbol{w}^T\boldsymbol{X}^T\boldsymbol{Xw} \\
&= \boldsymbol{y}^T \boldsymbol{y} \ -\ 2\boldsymbol{y}^T \boldsymbol{Xw}\ +\ \boldsymbol{w}^T\boldsymbol{X}^T\boldsymbol{Xw} \end{align*} $$

となります．

損失関数を最小化するｗを求める

さて、次はLを最小化するｗを求めるために、毎度おなじみの「(微分)＝0」を解きます。

Lをｗでベクトル微分して

$$\begin{align*}\frac{\partial L}{\partial w} &= \frac{\partial }{\partial w} ( \boldsymbol{y}^T \boldsymbol{y} \ -\ 2\boldsymbol{y}^T \boldsymbol{Xw}\ +\ \boldsymbol{w}^T\boldsymbol{X}^T\boldsymbol{Xw}) \\&= \boldsymbol{O} \ – (2\boldsymbol{y}^T \boldsymbol{X})^T + \{ \boldsymbol{X}^T\boldsymbol{X} +(\boldsymbol{X}^T\boldsymbol{X})^T\}\boldsymbol{w} \\&= -2\boldsymbol{X}^T \boldsymbol{y} + 2\boldsymbol{X}^T\boldsymbol{Xw} \end{align*}$$

よって、\(X^TX\)が正則なら

$$\begin{align*}\frac{\partial L}{\partial w} = 0 &\iff -2\boldsymbol{X}^T \boldsymbol{y} + 2\boldsymbol{X}^T\boldsymbol{Xw}=0 \\ &\iff \boldsymbol{X}^T\boldsymbol{Xw}=\boldsymbol{X}^T \boldsymbol{y} \\ &\iff \boldsymbol{w}=(\boldsymbol{X}^T \boldsymbol{X})^{-1} \boldsymbol{X}^T \boldsymbol{y} \end{align*}$$

となり、もとめる回帰係数ｗがわかりました。

＊なお \(X^TX\)が正則でないときは、\(X^TX\)に正則化項を加えて、正則になるように調整する方法が知られていますが、ここでは割愛します。気になる人はリッジ回帰やラッソ回帰について調べてみてください。

Pythonで実装

では、この結果 \(\boldsymbol{w}=(\boldsymbol{X}^T \boldsymbol{X})^{-1} \boldsymbol{X}^T \boldsymbol{y}\)をもちいて、Pythonで実装していきます。

全コードはこちら。

# ライブラリのインポート
import numpy as np
from sklearn.datasets import load_iris

# irisサンプルデータ
iris = load_iris()
X = iris.data[:,:-1]
y = iris.data[:,-1]

# 重回帰モデルの実装
class LinearRegression＿OLS():
    def __init__(self):
        self.coef_ = None
        
    def fit(self, X, y):
        X = np.insert(X, 0, 1, axis=1)
        self.coef_ = np.linalg.inv(X.T @ X) @ X.T @ y
    
    def predict(self, X):
        X = np.insert(X, 0, 1, axis=1)
        return X @ self.coef_

# 重回帰分析の実行
lr = LinearRegression_OLS()
lr.fit(X,y)
print(lr.predict(X))
print(lr.corf_)

必要なライブラリ・データの準備

# ライブラリのインポート
import numpy as np
from sklearn.datasets import load_iris

# irisサンプルデータ
iris = load_iris()
X = iris.data[:,:-1]
y = iris.data[:,-1]

まず、必要なライブラリをインポートします。今回はirisデータセットを使うために sklearn もインポートしました。

そしてXとyにはそれぞれ説明変数と目的変数を代入します。

今回の場合、Xには 150行3列 のデータが入っています。

fitメソッドの定義

# 重回帰モデルの実装
class LinearRegression＿OLS():
    def __init__(self):
        self.coef_ = None
        
    def fit(self, X, y):
        # Xの先頭列に x_0 = 1 を挿入してから計算
        X = np.insert(X, 0, 1, axis=1)
　　　　 #　w=(X^T X)^−1 X^T yを代入
        self.coef_ = np.linalg.inv(X.T @ X) @ X.T @ y

init関数では、coef_として回帰係数を定義しておきます。

そしてfitメソッドで、上で計算した \(\boldsymbol{w}=(\boldsymbol{X}^T \boldsymbol{X})^{-1} \boldsymbol{X}^T \boldsymbol{y}\) を変数coef_に代入します。逆行列はnumpy.linalg.invを使いましょう。行列の積の演算子は「@」です。

ここで一つ注意しなければならないことがあります。上の計算では、切片を入れるために便宜上の説明変数\(x_0\)を導入していたのでした。なので、何も処理せずにXをそのまま計算に使ってしまうとエラーが出るはずです。したがって、Xの先頭列に1の列を挿入してから計算しなければなりません。行列の挿入にはnumpyのinsertメソッドを使えばよいです。

predictメソッドの定義

    def predict(self, X):
        # Xの先頭列に x_0 = 1 を挿入してから計算
        X = np.insert(X, 0, 1, axis=1)
        # 予測値 y^ = Xw を返す
        return X @ self.coef_

回帰係数がcoef_に入っているので、Xにcoef_をかけてやれば予測値が求まります。

重回帰分析の実行

# 重回帰分析の実行
lr = LinearRegression_OLS()  # インスタンス化
lr.fit(X,y)  # 回帰係数の計算
print(lr.predict(X))  # 予測値の計算と出力
print(lr.corf_)  # 回帰係数の出力

実行の仕方は、scikit-learnのLinearRegressionのやつと全く同じです。

scikit-learnで実行した結果と比べてみても、ほぼ一致していることがわかると思います。