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

Python

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

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

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

  1. 最小二乗法で係数を算出 (線形代数を利用)
  2. その結果をもとにPythonで実装
  3. scikit-learnの結果と照らし合わせて確認

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

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

\(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*} $$

となります.

1行目→2行目:\(\hat{y}=Xw\)を代入

2行目→3行目:転置の公式 \((A+B)^T=A^T+B^T\)を適用

3行目→4行目: 転置の公式 \((AB)^T=B^TA^T\)を適用

4行目→5行目:分配法則で展開

5行目→6行目: \(\boldsymbol{y}^T \boldsymbol{Xw}\)と\(\boldsymbol{w}^T\boldsymbol{X}^T\boldsymbol{y}\)が転置の関係にあることを利用してまとめる

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

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

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

$$\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*}$$

1行目→2行目

第1項:スカラーなので微分するとゼロ(\(y^Ty\)はyのノルムの2乗)
第2項:微分公式 \(\frac{\partial}{\partial x}Ax = A^T\)を適用
第3項:微分公式 \(\frac{\partial}{\partial x}x^TAx = (A+A^T)x\)を適用

2行目→3行目:転置の公式 \((AB)^T=B^TA^T\)を適用してまとめる

よって、\(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*}$$

2行目→3行目: \(X^TX\)が正則なので、その逆行列を両辺の左側からかける

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

*なお \(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メソッドを使えばよいです。

numpy.insert(arr, obj, values, axis=None)

arr:もとの配列
obj:valuesを挿入したい位置[int, スライス,list]
values:挿入したい配列
axis:挿入する軸の方向

参考:https://numpy.org/doc/stable/reference/generated/numpy.insert.html

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で実行した結果と比べてみても、ほぼ一致していることがわかると思います。

ブログを読んでくれたあなたに質問です!

記事への意見・感想はコチラ