Skip to article frontmatterSkip to article content
Site not loading correctly?

This may be due to an incorrect BASE_URL configuration. See the MyST Documentation for reference.

線形回帰

モデル

線形回帰（linear regression）は、予測の目的変数yyと特徴量（説明変数）x1,x2,...,xdx_1, x_2, ..., x_dの間に次のような線形関係を仮定したモデルを置いて予測する手法。

y=β0+β1x1+⋯+βdxd+εy = \beta_0 + \beta_1 x_1 + \cdots + \beta_d x_d + \varepsilon
(1)

ここでβ1,β2,...,βd\beta_1, \beta_2, ..., \beta_dは回帰係数と呼ばれるパラメータで、モデル内で推定される。ε\varepsilonはデータ取得時の測定誤差などの偶然による誤差を表し、次の3つの条件を満たす。

  1. 期待値は0：E[ε]=0\text{E}[\varepsilon]=0

  2. 分散は一定: V[ε]=σ2\text{V}[\varepsilon]=\sigma^2

  3. 異なった誤差項は無相関: j≠ij \neq iならばCov(εi,εj)=E(εi,εj)=0Cov(\varepsilon_i, \varepsilon_j) = E(\varepsilon_i, \varepsilon_j) = 0

サンプルサイズがnnのデータセット{xi,yi}i=1n\{\boldsymbol{x}_i, y_i\}^n_{i=1}があるとして、目的変数をy=(y1,y2,...,yn)⊤\boldsymbol{y} = (y_1, y_2, ..., y_n)^\top、特徴量をX=(x1,x2,...,xn)⊤\boldsymbol{X}=(\boldsymbol{x}_1, \boldsymbol{x}_2, ..., \boldsymbol{x}_n)^\topとおくと、このモデルは

y=Xβ+ε\boldsymbol{y}=\boldsymbol{X} \boldsymbol{\beta}+\boldsymbol{\varepsilon}
(2)

と表記することができる。

パラメータの推定

一般的に線形回帰ではパラメータの推定に最小二乗法（least squares method）という方法が使われる。

これは誤差関数J(β)J( \boldsymbol{\beta})を実測値y\boldsymbol{y}と予測値y^=Xβ^\hat{\boldsymbol{y}} = \boldsymbol{X}\hat{\boldsymbol{\beta}}の二乗誤差の和（誤差二乗和 sum of squared error: SSE）

J(β)=∣∣y−Xβ^∣∣2=∑i=1n(yi−y^i)2=∑i=1nεi2=ε⊤εJ( \boldsymbol{\beta}) = ||\boldsymbol{y} -\boldsymbol{X}\hat{\boldsymbol{\beta}}||^2 = \sum^n_{i=1} (y_i - \hat{y}_i)^2 = \sum^n_{i=1} \varepsilon_i^2 = \boldsymbol{\varepsilon}^\top \boldsymbol{\varepsilon}
(3)

として定義し、この二乗誤差を最小にするパラメータ（最小二乗推定量 ordinary least square’s estimator: OLSE）

β^LS=arg minβ∑i=1n(yi−y^i)2\renewcommand{\argmin}{\mathop{\rm arg~min}\limits} \hat{\boldsymbol{\beta}}^{LS} = \argmin_{\boldsymbol{\beta}} \sum^n_{i=1} (y_i - \hat{y}_i)^2
(4)

を採用するという方法。

二乗誤差(yi−y^i)2=εi2(y_i - \hat{y}_i)^2 = \varepsilon_i^2はU字型になるため傾きがゼロになる点が最小値になる。そのため最小二乗法は解析的に解を求めることができる。

Source
import matplotlib.pyplot as plt
import numpy as np

def square(error):
    return error ** 2

errors = np.linspace(-1, 1, 100)
squared_errors = [square(error) for error in errors]

fig, ax = plt.subplots()
ax.plot(errors, squared_errors)
_ = ax.set(xlabel='error', ylabel='squared error')
<Figure size 432x288 with 1 Axes>

誤差二乗和は

ε⊤ε=(y−Xβ^)⊤(y−Xβ^)=y⊤y−y⊤Xβ−(Xβ)⊤y+(Xβ)⊤(Xβ)=y⊤y−2β⊤X⊤y+β⊤X⊤Xβ\begin{align} \boldsymbol{\varepsilon}^\top \boldsymbol{\varepsilon} &= (\boldsymbol{y} - \boldsymbol{X}\hat{\boldsymbol{\beta}})^\top (\boldsymbol{y} - \boldsymbol{X}\hat{\boldsymbol{\beta}})\\ &=\boldsymbol{y}^\top \boldsymbol{y} - \boldsymbol{y}^\top\boldsymbol{X}\boldsymbol{\beta} - (\boldsymbol{X}\boldsymbol{\beta})^\top\boldsymbol{y} + (\boldsymbol{X}\boldsymbol{\beta})^\top (\boldsymbol{X}\boldsymbol{\beta})\\ &= \boldsymbol{y}^\top \boldsymbol{y} - 2 \boldsymbol{\beta}^\top \boldsymbol{X}^\top \boldsymbol{y} + \boldsymbol{\beta}^\top \boldsymbol{X}^\top \boldsymbol{X} \boldsymbol{\beta}\\ \end{align}
(5)

であるから、二乗誤差の傾きがゼロになる点は

∂ε⊤ε∂β=−2X⊤y+2(X⊤X)β=0\frac{\partial \boldsymbol{\varepsilon}^\top \boldsymbol{\varepsilon}}{\partial \boldsymbol{\beta}} = -2\boldsymbol{X}^\top\boldsymbol{y} + 2(\boldsymbol{X}^\top\boldsymbol{X})\boldsymbol{\beta} =\boldsymbol{0}
(6)

と表すことができる。

これを整理して

2(X⊤X)β=2X⊤y2(\boldsymbol{X}^\top\boldsymbol{X}) \boldsymbol{\beta} = 2\boldsymbol{X}^\top\boldsymbol{y}
(7)

これの両辺を2で割ると（あるいは誤差関数の定義の際に1/21/2を掛けておくと）、正規方程式（normal equation）とよばれる次の式が得られる。

(X⊤X)β=X⊤y(\boldsymbol{X}^\top\boldsymbol{X})\boldsymbol{\beta} = \boldsymbol{X}^\top\boldsymbol{y}
(8)

これをβ\boldsymbol{\beta}について解けば

β=(X⊤X)−1X⊤y\boldsymbol{\beta} = (\boldsymbol{X}^\top\boldsymbol{X})^{-1} \boldsymbol{X}^\top\boldsymbol{y}
(9)

となり、最小二乗推定量β^LS\hat{\boldsymbol{\beta}}^{LS}が得られる。

実装

numpyでは、行列やベクトルの積は@という演算子で書くことができる。そのため、

import numpy as np
beta = np.linalg.inv(X.T @ X) @ X.T @ y

のように書けば上の式とおなじ演算を行うことができる。

データの準備

乱数を発生させて架空のデータを作る。

y=10+3x1+5x2+εx1∼Uniform(0,10)x2∼Normal(3,1)ε∼Normal(0,1)y = 10 + 3 x_1 + 5 x_2 + \varepsilon\\ x_1 \sim Uniform(0, 10)\\ x_2 \sim Normal(3, 1)\\ \varepsilon \sim Normal(0, 1)\\
(10)

ここでε\varepsilonは測定誤差などのランダムなノイズとする

import numpy as np
import pandas as pd
n = 100  # sample size

np.random.seed(0)
x0 = np.ones(shape=(n, ))
x1 = np.random.uniform(0, 10, size=n)
x2 = np.random.normal(3, 1, size=n)
noise = np.random.normal(size=n)

beta = [10, 3, 5]  # 真の回帰係数
y = beta[0] * x0 + beta[1] * x1 + beta[2] * x2 + noise

特徴量xxと目的変数yyの関係を散布図で描くと次の図のようになった。

Source
import matplotlib.pyplot as plt
import seaborn as sns

fig, axes = plt.subplots(ncols=2, figsize=(10, 3))
axes[0].scatter(x1, y)
axes[0].set(xlabel='x1', ylabel='y')
axes[1].scatter(x2, y)
axes[1].set(xlabel='x2', ylabel='y')
[Text(0.5, 0, 'x2'), Text(0, 0.5, 'y')]
<Figure size 720x216 with 2 Axes>

推定

これらのデータを使用して推定を行う。

X = np.array([x0, x1, x2]).T

# Xの冒頭5行は以下のようになっている
print(X[0:5])
[[1.         5.48813504 1.83485016]
 [1.         7.15189366 3.90082649]
 [1.         6.02763376 3.46566244]
 [1.         5.44883183 1.46375631]
 [1.         4.23654799 4.48825219]]

# 最小二乗法で推定
beta_ = np.linalg.inv(X.T @ X) @ X.T @ y

print(f"""
推定された回帰係数: {beta_.round(3)}
データ生成過程の係数: {beta}
""")

推定された回帰係数: [9.564 2.98  5.119]
データ生成過程の係数: [10, 3, 5]

真の値にそれなりに近い回帰係数が推定できた。

なお、scikit-learnに準拠したfit/predictのメソッドを持つ形でクラスとして定義するなら、以下のようになる（参考： sklearn準拠モデルの作り方 - Qiita）。

# scikit-learnに準拠した形で実装
from sklearn.base import BaseEstimator, RegressorMixin


class LinearRegression(BaseEstimator, RegressorMixin):

    def fit(self, X, y):
        self.coef_ = np.linalg.inv(X.T @ X) @ X.T @ y
        return self

    def predict(self, X):
        return X @ self.coef_

model = LinearRegression()
model.fit(X, y)
model.coef_
array([9.56372548, 2.97972446, 5.11931302])

予測してみる

root mean squared error (RMSE)

RMSE=1N∑i=1N(yi−y^i)2RMSE = \sqrt{ \frac{1}{N} \sum^N_{i=1} (y_i - \hat{y}_i)^2 }
(11)

を使って予測値を評価してみる。

# 予測値を算出
y_pred = model.predict(X)

# 予測誤差を評価
from sklearn.metrics import mean_squared_error
rmse = mean_squared_error(y, y_pred, squared=False)

print(f"RMSE: {rmse:.3f}")
RMSE: 0.972

予測値と実測値の散布図を描くと次のようになった。

Source
fig, ax = plt.subplots()
ax.plot(y_pred, y_pred, color='gray', alpha=.7)
ax.scatter(y_pred, y, alpha=.7)
_ = ax.set(xlabel='Predicted', ylabel='Actual')
[Text(0.5, 0, 'Predicted'), Text(0, 0.5, 'Actual')]
<Figure size 432x288 with 1 Axes>
機械学習
線形モデル
機械学習
Ridge回帰