線形回帰〜重回帰分析の理論とPythonのフルスクラッチ、scikit-learnによる実装〜

This article is available in: English

はじめに
多変数間の関係を捉える
最小二乗法を用いた回帰係数の推定
Pythonを用いて線形重回帰分析を行う
scikit-learnを用いて線形重回帰分析を行う

はじめに

　前回の記事では最小二乗法による線形回帰について述べました。

線形回帰〜最小二乗法の計算/解説とPythonのフルスクラッチ、scikit-learnよる実装〜

機械学習の中でも歴史ある線形回帰について、なかでも「最小二乗法」についての理論とpython, scikit-learnによる実装を紹介します

そこでは、バネに吊るした重りの質量を説明変数、バネの長さを目的変数として回帰直線を求めました。つまり、説明変数が１つのみの場合でしたが、では説明変数が複数ある場合にはどのように回帰モデルを構成すれば良いのでしょうか？

具体例として、アパートの家賃を挙げてみます。表１のように、
・面積　・築年数　・駅徒歩分の３つの説明変数、
・家賃を目的変数
とするデータです。面積、築年数、駅徒歩分のデータを用いて家賃の予測モデルを構築するには、どういった定式化を行えば良いでしょう？

▲表1: 3つの説明変数を用いて、目的変数を予測するモデルを構築したい

　この記事では、複数の説明変数と目的変数の関係をモデル化する重回帰分析と、回帰係数を求める正規方程式の導出を行います。

また、Pythonによる重回帰分析の実装も行います。コードは以下で試すことができます。

Google Colaboratory

\begin{align*}
\newcommand{\mat}[1]{\begin{pmatrix} #1 \end{pmatrix}}
\newcommand{\f}[2]{\frac{#1}{#2}}
\newcommand{\pd}[2]{\frac{\partial #1}{\partial #2}}
\newcommand{\d}[2]{\frac{{\rm d}#1}{{\rm d}#2}}
\newcommand{\T}{\mathsf{T}}
\newcommand{$}{\left(}
\newcommand{$}{\right)}
\newcommand{\{}{\left\{}
\newcommand{\}}{\right\}}
\newcommand{\[}{\left[}
\newcommand{\]}{\right]}
\newcommand{\dis}{\displaystyle}
\newcommand{\eq}[1]{{\rm Eq}(\ref{#1})}
\newcommand{\n}{\notag\\}
\newcommand{\t}{\ \ \ \ }
\newcommand{\tt}{\t\t\t\t}
\newcommand{\argmax}{\mathop{\rm arg\, max}\limits}
\newcommand{\argmin}{\mathop{\rm arg\, min}\limits}
\def\l<#1>{\left\langle #1 \right\rangle}
\def\us#1_#2{\underset{#2}{#1}}
\def\os#1^#2{\overset{#2}{#1}}
\newcommand{\case}[1]{\{ \begin{array}{ll} #1 \end{array} \right.}
\newcommand{\s}[1]{{\scriptstyle #1}}
\definecolor{myblack}{rgb}{0.27,0.27,0.27}
\definecolor{myred}{rgb}{0.78,0.24,0.18}
\definecolor{myblue}{rgb}{0.0,0.443,0.737}
\definecolor{myyellow}{rgb}{1.0,0.82,0.165}
\definecolor{mygreen}{rgb}{0.24,0.47,0.44}
\newcommand{\c}[2]{\textcolor{#1}{#2}}
\newcommand{\ub}[2]{\underbrace{#1}_{#2}}
\end{align*}

多変数間の関係を捉える

　上記では例として説明変数を3つとしましたが、もう少し一般的に $p$ 個の説明変数 $x_1, x_2, \dots, x_p$ と目的関数 $y$ について考えます。

以下の表2は表1を一般化したもので、$n$ 個データを測定した結果を示します。例えば $i$ 番目のデータ点では説明変数の組 $x_1^{(i)}, x_2^{(i)}, \dots, x_p^{(i)}$ に対して目的変数 $y^{(i)}$ が関係しています。

▲表2: 説明変数が $p$ 個のデータを $n$ 個観測した際の対応表。列が説明変数の種類を表し、行がデータのindexを表す

目的は、目的変数 $y$ と説明変数 $x_1, x_2, \dots, x_p$ を結びつけるモデルを構築することにあります。ここでは、データの傾向や、先験的な知見からこの現象の構造は線形であると仮定します。

すなわち、目的変数と説明変数の関係性は、線形式

\begin{align*}
y = \theta_0 + \theta_1 x_1 + \theta_1 x_2 + \cdots + \theta_p x_p
\end{align*}

で表現できるとします。ここで、$\theta_0, \dots, \theta_p$ は回帰係数です。

実際には誤差 $\varepsilon$ の影響でデータは回帰式からズレてしまっていると考えると、各データ点は次の関係を満たします。（データ番号の添字を右上に付けることにします）

\begin{align*}
y^{(1)} &= \theta_0 + \theta_1 x_1^{(1)} + \theta_1 x_2^{(1)} + \cdots + \theta_p x_p^{(1)} + \varepsilon ^{(1)}, \n
y^{(2)} &= \theta_0 + \theta_1 x_1^{(2)} + \theta_1 x_2^{(2)} + \cdots + \theta_p x_p^{(2)} + \varepsilon ^{(2)}, \n
& \ \ \vdots \n
y^{(i)} &= \theta_0 + \theta_1 x_1^{(i)} + \theta_1 x_2^{(i)} + \cdots + \theta_p x_p^{(i)} + \varepsilon ^{(i)}, \n
&\ \ \vdots \n
y^{(n)} &= \theta_0 + \theta_1 x_1^{(n)} + \theta_1 x_2^{(n)} + \cdots + \theta_p x_p^{(n)} + \varepsilon ^{(n)}.
\end{align*}

このように2つ以上の説明変数と目的変数との関係を捉えるモデルは、線形重回帰モデルと呼ばれます。

　さて、上式は行列を用いて表現すると、次のように表すことができます。

\begin{align*}
\underbrace{\mat{y^{(1)} \\ y^{(2)} \\ \vdots \\ y^{(i)} \\ \vdots \\y^{(n)}}}_{\equiv \bm{y}}
=
\underbrace{\mat{
1 & x_1^{(1)} & x_2^{(1)} & \cdots & x_p^{(1)} \\
1 & x_1^{(2)} & x_2^{(2)} & \cdots & x_p^{(2)} \\
\vdots & \vdots & \vdots & \ddots & \vdots \\
1 & x_1^{(i)} & x_2^{(i)} & \cdots & x_p^{(i)} \\
\vdots & \vdots & \vdots & \ddots & \vdots \\
1 & x_1^{(n)} & x_2^{(n)} & \cdots & x_p^{(n)}
}}_{\equiv X}
\underbrace{\mat{\theta_0 \\ \theta_1 \\ \theta_2 \\ \vdots \\ \theta_p}}_{\equiv \bm{\theta}}
+
\underbrace{\mat{\varepsilon^{(1)} \\ \varepsilon^{(2)} \\ \vdots \\ \varepsilon^{(i)} \\ \vdots \\\varepsilon^{(n)}}}_{\equiv \bm{\varepsilon}}.
\end{align*}

さらに、この式をベクトルと行列を用いて表すと、

\begin{align}
\bm{y} = X \bm{\theta} + \bm{\varepsilon}
\label{eq:reg}
\end{align}

となることがわかります。

では次節でこの回帰係数ベクトル $\bm{\theta}$ を決定する方法について述べていきます。

最小二乗法を用いた回帰係数の推定

　説明変数が１変数の場合と同様に、最小二乗法を用いて回帰係数ベクトルを推定します。

つまり、各データでの誤差項の２乗和 $J \equiv (\varepsilon^{(1)})^2 + (\varepsilon^{(2)})^2 + \cdots + (\varepsilon^{(n)})^2$ を最小とするように回帰係数ベクトル $\bm{\theta}$ の推定値 $\hat{\bm{\theta}}$ を求める方法です。

$\eq{eq:reg}$ より、$\bm{\varepsilon} = \bm{y}\, – X \bm{\theta}$ となるので、線形重回帰での最小二乗法を数式で表現すると、以下となります。

\begin{align*}
\hat{\bm{\theta}} &= \argmin_{\bm{\theta}} J(\bm{\theta}) \n
&= \argmin_{\bm{\theta}} \sum_{i=1}^n (\varepsilon^{(i)})^2 \n
&= \argmin_{\bm{\theta}}\, \bm{\varepsilon}^\T \bm{\varepsilon} \n
&= \argmin_{\bm{\theta}}\, (\bm{y}\, – X \bm{\theta})^\T (\bm{y}\, – X \bm{\theta}).
\end{align*}

したがって、前回と同様に $J(\bm{\theta})$ を $\bm{\theta}$ で偏微分し、$\bm{0}$ と置くことで推定値を求めます。

\begin{align*}
J(\bm{\theta}) &= (\bm{y}\, – X \bm{\theta})^\T (\bm{y}\, – X \bm{\theta}) \n
&= (\bm{y}^\T – \bm{\theta}^\T X^\T) (\bm{y}\, – X \bm{\theta}) \t \t (\because (AB)^\T = B^\T A^\T) \n
&= \bm{y}^\T \bm{y}\ – \bm{y}^\T X \bm{\theta}\ – \underbrace{\bm{\theta}^\T X^\T \bm{y}}_{= \bm{y}^\T X \bm{\theta}} + \bm{\theta}^\T X^\T X \bm{\theta} \n
&= \bm{y}^\T \bm{y}\ – 2 \bm{y}^\T X \bm{\theta} + \bm{\theta}^\T X^\T X \bm{\theta}
\end{align*}

ですから、ベクトル微分の公式を使って極値を求めます。

\begin{align*}
\pd{J(\bm{\theta})}{\bm{\theta}} &= \pd{}{\bm{\theta}} $ \bm{y}^\T \bm{y}\ – 2 \bm{y}^\T X \bm{\theta} + \bm{\theta}^\T X^\T X \bm{\theta} $ \n
&= \underbrace{\pd{}{\bm{\theta}} (\bm{y}^\T \bm{y})}_{=0} \ – 2 \pd{}{\bm{\theta}} (\bm{y}^\T X \bm{\theta}) + \pd{}{\bm{\theta}} (\bm{\theta}^\T X^\T X \bm{\theta}) \n
&= -2 X^\T \bm{y} + 2X^\T X \bm{\theta} = \bm{0}.
\end{align*}

ここで、最後の行の変形ではベクトル微分の公式を用いました。

「ベクトルで微分」公式導出

機械学習の理論を勉強すると「スカラーをベクトルで微分」する操作をよく目にします。この記事では「スカラーをベクトルで微分」の公式を導出します。

したがって、正規方程式と呼ばれる次式を得ます。

\begin{align*}
X^\T X \bm{\theta} = X^\T \bm{y}
\end{align*}

以上より、もし行列 $X^\T X$ の逆行列が存在すれば、最小二乗推定値は下記であたえられます。

\begin{align*}
\hat{\bm{\theta}} = (X^\T X)^{-1} X^\T \bm{y}\t (☆)
\end{align*}

Pythonを用いて線形重回帰分析を行う

　それでは、Pythonをもちいて線形重回帰分析をおこないましょう。式(☆)を計算すれば回帰係数ベクトル $\hat{\bm{\theta}}$ が求まります。

例として用いるデータは、可視化の都合から、説明変数が2個とします。

import numpy as np
import pandas as pd

area = [17.00, 18.00, 21.00, 24.00, 18.90, 20.66, 23.51, 25.21, 24.94, 30.22]  # 説明変数1
age = [31, 34, 36, 31, 28, 22, 19, 12, 4, 0]  # 説明変数2
rent = [35000, 36000, 39000, 42000, 46000, 50000, 59000, 68000, 73000, 86000]  # 目的変数

df = pd.DataFrame({
    'area': area,
    'age': age,
    'rent': rent})

まずはplotlyを用いて、対象とする3次元データの可視化を行います。

import plotly.express as px
import plotly.graph_objects as go

# データを可視化してみる
fig = px.scatter_3d(df, x='area', y='age', z='rent', color='rent', opacity=0.7)
fig.show()

このデータに合う回帰平面を求めます。

では、式 (☆) を計算するために行列 $X$ 、ベクトル $\bm{y}$ を用意します。

# 行列X, ベクトルyの準備

X = df[['area', 'age']].values
ones = np.ones(len(X)).reshape(-1, 1)
X = np.hstack((ones, X))
print('行列X')
print(X)

y = df['rent'].values
print('ベクトルy')
print(y)

出力
行列X
[[ 1.   17.   31.  ]
 [ 1.   18.   34.  ]
 [ 1.   21.   36.  ]
 [ 1.   24.   31.  ]
 [ 1.   18.9  28.  ]
 [ 1.   20.66 22.  ]
 [ 1.   23.51 19.  ]
 [ 1.   25.21 12.  ]
 [ 1.   24.94  4.  ]
 [ 1.   30.22  0.  ]]
ベクトルy
[35000 36000 39000 42000 46000 50000 59000 68000 73000 86000]

式 (☆)を計算して回帰係数ベクトル $\hat{\bm{\theta}}$ を求めます。

def multiple_regression(X, y):
    """回帰係数ベクトルを計算する"""
    theta = np.linalg.pinv(X.T @ X) @ X.T @ y
    return theta

theta = multiple_regression(X, y)
theta_0 = theta[0]
theta_1 = theta[1]
theta_2 = theta[2]

print(f'theta_0: {theta_0}')
print(f'theta_1: {theta_1}')
print(f'theta_2: {theta_2}')

出力
theta_0: 50300.79618461159
theta_1: 1154.3386202815495
theta_2: -1045.7759592710272

回帰係数が求まったので、回帰平面をプロットしてみましょう。

# 回帰平面の可視化
mesh_size = 1
margin = 0.1
x1_min, x1_max = df['area'].min()-margin, df['area'].max()+margin
x2_min, x2_max = df['age'].min()-margin, df['age'].max()+margin
x1_range = np.arange(x1_min, x1_max, mesh_size)
x2_range = np.arange(x2_min, x2_max, mesh_size)
xx1, xx2 = np.meshgrid(x1_range, x2_range)

y = (theta_0 + theta_1 * xx1 + theta_2 * xx2)  # 回帰平面

fig = px.scatter_3d(df, x='area', y='age', z='rent', opacity=0.8)  # データ点のプロット
fig.update_traces(marker=dict(
    color=colors[0],
    size=8,
    line=dict(width=2,color='white')))

fig.add_traces(go.Surface(x=x1_range, y=x2_range, z=y, opacity=0.7))  # 平面のプロット
fig.show()

このように、正規方程式を解くことで回帰平面を得ることができました。

scikit-learnを用いて線形重回帰分析を行う

　scikit-learnでも同様に重回帰分析を行なうことができます。

from sklearn.linear_model import LinearRegression

area = [17.00, 18.00, 21.00, 24.00, 18.90, 20.66, 23.51, 25.21, 24.94, 30.22]  # 説明変数1
age = [31, 34, 36, 31, 28, 22, 19, 12, 4, 0]  # 説明変数2
rent = [35000, 36000, 39000, 42000, 46000, 50000, 59000, 68000, 73000, 86000]  # 目的変数
df = pd.DataFrame({
    'area': area,
    'age': age,
    'rent': rent})

X = df[['area', 'age']].values
y = df['rent'].values

# 回帰係数を求める
model = LinearRegression()
model.fit(X, y)

theta_0 = model.intercept_
theta_1 = model.coef_[0]
theta_2 = model.coef_[1]

print(f'theta_0: {theta_0}')
print(f'theta_1: {theta_1}')
print(f'theta_2: {theta_2}')

出力
theta_0: 50300.79618451819
theta_1: 1154.3386202847946
theta_2: -1045.7759592701218

# 回帰平面の可視化
mesh_size = 1
margin = 0.1
x1_min, x1_max = df['area'].min()-margin, df['area'].max()+margin
x2_min, x2_max = df['age'].min()-margin, df['age'].max()+margin
x1_range = np.arange(x1_min, x1_max, mesh_size)
x2_range = np.arange(x2_min, x2_max, mesh_size)
xx1, xx2 = np.meshgrid(x1_range, x2_range)

y = (theta_0 + theta_1 * xx1 + theta_2 * xx2)  # 回帰平面

fig = px.scatter_3d(df, x='area', y='age', z='rent', opacity=0.8)  # データ点のプロット
fig.update_traces(marker=dict(
    color=colors[0],
    size=8,
    line=dict(width=2,color='white')))

fig.add_traces(go.Surface(x=x1_range, y=x2_range, z=y, opacity=0.7))  # 平面のプロット
fig.show()

scikit-learnでも同様に回帰平面を得ることができました。

リンク