2. 重回帰#

!pip install japanize-matplotlib
Requirement already satisfied: japanize-matplotlib in /home/okazaki/.local/lib/python3.8/site-packages (1.1.3)
Requirement already satisfied: matplotlib in /usr/local/lib/python3.8/dist-packages (from japanize-matplotlib) (3.5.1)
Requirement already satisfied: numpy>=1.17 in /usr/local/lib/python3.8/dist-packages (from matplotlib->japanize-matplotlib) (1.22.2)
Requirement already satisfied: cycler>=0.10 in /usr/local/lib/python3.8/dist-packages (from matplotlib->japanize-matplotlib) (0.11.0)
Requirement already satisfied: fonttools>=4.22.0 in /usr/local/lib/python3.8/dist-packages (from matplotlib->japanize-matplotlib) (4.29.1)
Requirement already satisfied: python-dateutil>=2.7 in /usr/local/lib/python3.8/dist-packages (from matplotlib->japanize-matplotlib) (2.8.2)
Requirement already satisfied: pyparsing>=2.2.1 in /usr/local/lib/python3.8/dist-packages (from matplotlib->japanize-matplotlib) (3.0.7)
Requirement already satisfied: kiwisolver>=1.0.1 in /usr/local/lib/python3.8/dist-packages (from matplotlib->japanize-matplotlib) (1.3.2)
Requirement already satisfied: pillow>=6.2.0 in /usr/local/lib/python3.8/dist-packages (from matplotlib->japanize-matplotlib) (9.0.1)
Requirement already satisfied: packaging>=20.0 in /usr/local/lib/python3.8/dist-packages (from matplotlib->japanize-matplotlib) (21.3)
Requirement already satisfied: six>=1.5 in /usr/lib/python3/dist-packages (from python-dateutil>=2.7->matplotlib->japanize-matplotlib) (1.14.0)
import numpy as np
import matplotlib.pyplot as plt
import japanize_matplotlib

2.1. 重回帰#

これまで、以下のデータ\(\mathcal{D}_s\)を例として単回帰を説明した。

(2.1)#\[\begin{align} \mathcal{D}_s = \left\{(x_i, y_i)\right\}_{i=1}^{4} = \left\{(1, 3), (3, 6), (6, 5), (8, 7)\right\} \end{align}\]

また、求めた回帰直線をデータ点とともにプロットしたものを以下に示す。

a, b = 0.431, 3.310 
D = np.array([[1, 3], [3, 6], [6, 5], [8, 7]])

fig, ax = plt.subplots(dpi=100, figsize=(6, 6))

ax.scatter(D[:,0], D[:,1], marker='o')

x = np.array([0, 10])
ax.plot(x, a * x + b, 'tab:red', ls='-', label=f'${a}x + {b}$')

for i, row in enumerate(D):
    x, y = row
    y_hat = a * x + b
    plt.vlines([x], min(y, y_hat), max(y, y_hat), "tab:blue", linestyles='dashed')

ax.set_xlabel('$x$')
ax.set_ylabel('$y$')
ax.set_xlim(0, 10)
ax.set_ylim(0, 10)
ax.set_aspect('equal')
ax.grid()
plt.legend(loc='upper left')
plt.show()
../_images/02mra_4_0.png

残差の平均二乗和が最小になるような直線を引いたはずであるが、残差を完全に無くすことはできなかった(決定係数を求めてみると\(0.616\)であり、あまり高くない)。プロットされている点を眺めると、3次関数を当てはめる方が向いているように思われる。そこで、1次関数よりも表現力が高いモデルとして、\(d\)次の多項式の一般形、

(2.2)#\[\begin{align} \hat{y} = w_0 + w_1 x + w_2 x^2 + \dots + w_d x^d \end{align}\]

にフィッティングすることを考えたい。ここで、\(w_0, w_1, w_2, \dots, w_d \in \mathbb{R}\)はモデルのパラメータである。

なお、\(\hat{y}\)は重み\(w_0, w_1, w_2, \dots, w_d\)による説明変数\(1, x, x^2, \dots x^d\)の線形重み付き和と考えることができる。このモデルは多項式の一般形ではあるが、\(\hat{y}\)\(w_i\)\(x^i\)の線形結合で表されているため、線形モデルと見なすことができる。ただし、説明変数が2個以上であるので、単回帰を重回帰に拡張する必要がある。

一般的に、\(d\)個の説明変数\(x_1, x_2, \dots, x_d\)に対して、線形モデル、すなわち\(d\)個の重み\(w_1, w_2, \dots, w_d\)の線形結合で目的変数の予測値\(\hat{y}\)を計算することを線形回帰と呼ぶ。

(2.3)#\[\begin{align} \hat{y} = w_1 x_1 + w_2 x_2 + \dots + w_d x_d \end{align}\]

ゆえに、\((d-1)\)次の多項式へのフィッティングは、\(d\)個の説明変数\(x_1 = 1\), \(x_2 = x\), \(x_3 = x^2\), \(\dots\), \(x_d = x^{d-1}\)による線形回帰と見なすことができる。下図は、重回帰でデータ\(\mathcal{D}_s\)に3次関数を当てはめた結果である。本章の内容をマスターすると、\(4\)個の説明変数\(1, x, x^2, x^3\)に対して、その重み(\(-1.23, 5.41, -1.27, 0.09\))を求めることができる。

a, b = 0.431, 3.310 
D = np.array([[1, 3], [3, 6], [6, 5], [8, 7]])

fig, ax = plt.subplots(dpi=100, figsize=(6, 6))

ax.scatter(D[:,0], D[:,1], marker='o')

x = np.linspace(0, 10, 500)
ax.plot(x, 0.09 * x ** 3 - 1.27 * x ** 2 + 5.41 * x - 1.23, 'tab:red', ls='-', label='$0.09 x^3 - 1.27 x^2 + 5.41 x - 1.23$')
ax.set_xlabel('$x$')
ax.set_ylabel('$y$')
ax.set_xlim(0, 10)
ax.set_ylim(0, 10)
ax.set_aspect('equal')
ax.grid()
plt.legend(loc='upper left')
plt.show()
../_images/02mra_6_0.png

回帰分析において、説明変数を何個用意するか、それぞれの説明変数をどのように定義するかは任意である。たとえば、東京におけるアイスクリームへの支出額を予測するために、最高気温だけでなく、予測したい月の数字やアイスクリームの平均価格も有力な説明変数になり得るかもしれない。また、最高気温の平方根や、最高気温と月の積など、元の説明変数に非線形変換を適用した説明変数を導入してもよい。分析者が非線形変換の方法を明示的に設計する必要があるが、説明変数の組み合わせを考慮した非線形な回帰分析も、線形回帰の枠組みで行うことができる。

non-linear

2.2. 表記法の定義#

重回帰の目的関数を一般的に記述するため、ベクトルによる表記を導入する。\(d\)個の説明変数\(x_1, x_2, \dots, x_d\)をベクトル\(\pmb{x} \in \mathbb{R}^d\)でまとめて表現する。

(2.4)#\[\begin{align} \pmb{x} = \begin{pmatrix} x_1 \\ x_2 \\ \vdots \\ x_d \end{pmatrix} \end{align}\]

また、\(d\)個のパラメータ(重み)\(w_1, w_2, \dots, w_d\)をベクトル\(\pmb{w} \in \mathbb{R}^d\)でまとめて表現する。

(2.5)#\[\begin{align} \pmb{w} = \begin{pmatrix} w_1 \\ w_2 \\ \vdots \\ w_d \end{pmatrix} \end{align}\]

すると、重回帰モデルによる目的変数の予測値\(\hat{y}\)は、ベクトル\(\pmb{x}\)\(\pmb{w}\)の内積として記述できる。

(2.6)#\[\begin{align} \hat{y} = \pmb{x}^\top\pmb{w} = w_1 x_1 + w_2 x_2 + \dots + w_d x_d \end{align}\]

単回帰のときと同様に、説明変数と目的変数の組を事例として表現する。単回帰の場合と異なるのは、説明変数がスカラー\(x\)からベクトル\(\pmb{x}\)に拡張されたことである。 \(1\)番目の事例を\((\pmb{x}_1, y_1)\)\(2\)番目の事例を\((\pmb{x}_2, y_2)\)\(i\)番目の事例を\((\pmb{x}_i, y_i)\)と表すことにすると、\(N\)個の事例からなるデータ\(\mathcal{D}\)は次のように表される。

(2.7)#\[\begin{align} \mathcal{D} = \left\{(\pmb{x}_1, y_1), (\pmb{x}_2, y_2), \dots, (\pmb{x}_N, y_N)\right\} = \left\{(\pmb{x}_i, y_i)\right\}_{i=1}^{N} \end{align}\]

さらに、\(N\)個の事例や目的変数を行列やベクトルで記述する表記法も導入しておく。説明変数のベクトル\(\pmb{x}_i\)を行ベクトルとして縦に\(N\)個並べた行列\(\pmb{X} \in \mathbb{R}^{N \times d}\)を定義する。この行列\(\pmb{X}\)計画行列(design matrix)と呼ばれ、それぞれの事例の\(d\)個の説明変数を水平方向に\(d\)個並べ、\(N\)個の事例を垂直方向に\(N\)個並べたものである。

(2.8)#\[\begin{align} \pmb{X} = \begin{pmatrix} \pmb{x}_1^\top \\ \pmb{x}_2^\top \\ \vdots \\ \pmb{x}_N^\top \\ \end{pmatrix} = \begin{pmatrix} X_{1,1} & X_{1,2} & \dots & X_{1,d} \\ X_{2,1} & X_{2,2} & \dots & X_{2,d} \\ \vdots & \vdots & \ddots & \vdots \\ X_{N,1} & X_{N,2} & \dots & X_{N,d} \end{pmatrix} \end{align}\]

ここで、\(X_{i,j}\)は行列\(\pmb{X}\)\(i\)\(j\)列の要素を表す。

また、目的変数\(y_i\)を縦に\(N\)個並べたベクトル\(\pmb{y} \in \mathbb{R}^N\)を定義する。

(2.9)#\[\begin{align} \pmb{y} = \begin{pmatrix} y_1 \\ y_2 \\ \vdots \\ y_N \end{pmatrix} \end{align}\]

目的変数の予測値\(\hat{y}_i\)を縦に\(N\)個並べ、ベクトル\(\hat{\pmb{y}} \in \mathbb{R}^N\)と書くことにすると、その予測値は計画行列とパラメータベクトルの積で求めることができる。

(2.10)#\[\begin{align} \hat{\pmb{y}} = \begin{pmatrix} \hat{y}_1 \\ \hat{y}_2 \\ \vdots \\ \hat{y}_N \end{pmatrix} = \begin{pmatrix} \pmb{x}_1^\top \pmb{w} \\ \pmb{x}_2^\top \pmb{w} \\ \vdots \\ \pmb{x}_N^\top \pmb{w} \end{pmatrix} = \pmb{X} \pmb{w} \end{align}\]

2.2.1. 単回帰を表現する#

これまでに定義した表記法で単回帰を記述できることを確認しておこう。\(d=2\)として説明変数とパラメータベクトルを以下のように定義する。

(2.11)#\[\begin{split} \begin{align} \pmb{x} = \begin{pmatrix} 1 \\ x \end{pmatrix}, \pmb{w} = \begin{pmatrix} b \\ a \end{pmatrix} \end{align} \end{split}\]

目的変数の予測値\(\hat{y}\)を計算する式を展開すると、単回帰モデルが得られる。

(2.12)#\[\begin{align} \hat{y} = \pmb{x}^\top\pmb{w} = b + xa = ax + b \end{align}\]

冒頭のグラフの例(\(N=4\))の学習データは次のように表現される。

\[\begin{align*} \mathcal{D}_s = \left\{(\pmb{x}_i, y_i)\right\}_{i=1}^{4} = \left\{\left(\begin{pmatrix}1 & 1\end{pmatrix}^\top, 3\right), \left(\begin{pmatrix}1 & 3\end{pmatrix}^\top, 6\right), \left(\begin{pmatrix}1 & 6\end{pmatrix}^\top, 5\right), \left(\begin{pmatrix}1 & 8\end{pmatrix}^\top, 7\right)\right\} \end{align*}\]

計画行列\(\pmb{X} \in \mathbb{R}^{4 \times 2}\)と目的変数ベクトル\(\pmb{y} \in \mathbb{R}^4\)は、それぞれ、

(2.13)#\[\begin{align} \pmb{X} = \begin{pmatrix} 1 & 1 \\ 1 & 3 \\ 1 & 6 \\ 1 & 8 \end{pmatrix}, \pmb{y} = \begin{pmatrix} 3 \\ 6 \\ 5 \\ 7 \end{pmatrix} \end{align}\]

目的変数の予測ベクトル\(\hat{\pmb{y}}\)は行列\(W\)とベクトル\(\pmb{w}\)の積で求めることができる。

(2.14)#\[\begin{align} \hat{\pmb{y}} = \pmb{X} \pmb{w} = \begin{pmatrix} 1 & 1 \\ 1 & 3 \\ 1 & 6 \\ 1 & 8 \end{pmatrix} \begin{pmatrix} b \\ a \end{pmatrix} = \begin{pmatrix} a + b \\ 3a + b \\ 6a + b \\ 8a + b \end{pmatrix} \end{align}\]

2.3. 目的関数#

単回帰のときと同様に、ある事例の目的変数の真の値\(y_i\)と推定値\(\hat{y}_i\)の差を残差\(\epsilon_i\)とする。

(2.15)#\[\begin{align} \epsilon_i = y_i - \hat{y}_i = y_i - \pmb{x}_i^\top\pmb{w} \end{align}\]

学習データにおける平均二乗残差は、

(2.16)#\[\begin{align} \hat{L} = \frac{1}{N} \sum_{i=1}^{N} \epsilon_i^2 = \frac{1}{N} \sum_{i=1}^{N} (y_i - \hat{y}_i)^2 = \frac{1}{N} \sum_{i=1}^{N} (y_i - \pmb{x}_i^\top\pmb{w})^2 \end{align}\]

重回帰モデルの学習において平均二乗残差を最小化するときは、学習データ\(\mathcal{D}\)や学習事例数\(N\)は事前に与えられる定数であり、パラメータ\(\pmb{w}\)は変数である。したがって、\(N\)は最小化の解に影響を与えないので(二乗残差の平均を最小化することと二乗残差の和を最小化することは同じ)、今後は目的関数に含めないことにする。すると、重回帰で最小化したい目的関数は、

\[ \begin{align} \hat{L}_{\mathcal{D}}(\pmb{w}) = \sum_{i=1}^{N} (y_i - \pmb{x}_i^\top\pmb{w})^2 \end{align} \]

なお、全事例に関する残差をベクトル、

(2.17)#\[\begin{align} \pmb{\epsilon} = \begin{pmatrix} \epsilon_1 \\ \epsilon_2 \\ \vdots \\ \epsilon_N \end{pmatrix} \end{align}\]

で表すことにすると、

(2.18)#\[ \begin{align} \pmb{\epsilon} = \pmb{y} - \hat{\pmb{y}} = \pmb{y} - \pmb{X}\pmb{w} \end{align} \]

である。

以上まとめると、重回帰で最小化したい目的関数は、

重回帰の目的関数(二乗残差和)

(2.19)#\[ \begin{align} \hat{L}_{\mathcal{D}}(\pmb{w}) = \lVert \pmb{\epsilon} \rVert^2 = \lVert \pmb{y} - \pmb{X}\pmb{w} \rVert^2 = \sum_{i=1}^{N} (y_i - \pmb{x}_i^\top\pmb{w})^2 \end{align} \]

と表現できる。

2.4. 目的関数の最小化#

目的関数\(\hat{L}_{\mathcal{D}}(\pmb{w})\)を最小にするパラメータのベクトル\(\pmb{w}\)を求めるため、式(2.19)をあるパラメータ\(w_j\)で偏微分する(\(j \in \{1, 2, \dots, d\}\)はパラメータのインデックスである)。

(2.20)#\[\begin{split} \begin{align} \frac{\partial \hat{L}_{\mathcal{D}}(\pmb{w})}{\partial w_j} &= \sum_{i=1}^N 2 \cdot \underbrace{(y_i - \pmb{x}_i^\top\pmb{w})}_{=\epsilon_i} \cdot \{-\underbrace{(\pmb{x}_i)_j}_{=X_{i,j}}\} \\ &= -2 \sum_{i=1}^N X_{i,j} \epsilon_i \\ &= -2 \begin{pmatrix} X_{1,j} & X_{2,j} & \dots & X_{N,j} \end{pmatrix} \begin{pmatrix} \epsilon_1 \\ \epsilon_2 \\ \vdots \\ \epsilon_N \end{pmatrix} \\ &= -2 (\pmb{X}^\top)_j \pmb{\epsilon} \\ &= -2 (\pmb{X}^\top \pmb{\epsilon})_j \end{align} \end{split}\]

なお、最後から2番目の式変形は、計画行列\(\pmb{X}\)の転置\(\pmb{X}^\top\)

(2.21)#\[\begin{align} \pmb{X}^\top = \begin{pmatrix} X_{1,1} & X_{2,1} & \dots & X_{N,1} \\ X_{1,2} & X_{2,2} & \dots & X_{N,2} \\ \vdots & \vdots & \ddots & \vdots \\ X_{1,j} & X_{2,j} & \dots & X_{N,j} \\ \vdots & \vdots & \ddots & \vdots \\ X_{1,d} & X_{2,d} & \dots & X_{N,d} \end{pmatrix} \end{align}\]

\(j\)行目の行ベクトルが\(\begin{pmatrix} X_{1,j} & X_{2,j} & \dots & X_{N,j}\end{pmatrix}\)であることを利用した。最終行の式変形では、以下の行列ベクトル積の\(j\)行目の計算結果に基づく。

(2.22)#\[\begin{align} \pmb{X}^\top \pmb{\epsilon} = \begin{pmatrix} (\pmb{X}^\top)_{1} \pmb{\epsilon} \\ (\pmb{X}^\top)_{2} \pmb{\epsilon} \\ \vdots \\ (\pmb{X}^\top)_{j} \pmb{\epsilon} \\ \vdots \\ (\pmb{X}^\top)_{d} \pmb{\epsilon} \\ \end{pmatrix} \end{align}\]

(2.20)より、\(\frac{\partial \hat{L}_{\mathcal{D}}(\pmb{w})}{\partial w_j}\)\(-2 (\pmb{X}^\top \pmb{\epsilon})\)\(j\)行目の要素に対応する。全ての\(j \in \{1, 2, \dots, d\}\)について偏微分したものを\(d\)次元ベクトルとしてまとめ、\(\nabla \hat{L}_{\mathcal{D}}(\pmb{w})\)と書くことにすると、

\[\begin{split} \nabla \hat{L}_{\mathcal{D}}(\pmb{w}) = \begin{pmatrix} \frac{\partial }{\partial w_1} \\ \frac{\partial }{\partial w_2} \\ \vdots \\ \frac{\partial }{\partial w_d} \\ \end{pmatrix} \hat{L}_{\mathcal{D}}(\pmb{w}) = -2 \pmb{X}^\top \pmb{\epsilon} \end{split}\]

さらに、式(2.18)より\(\pmb{\epsilon} = \pmb{y} - \pmb{X}\pmb{w}\)であるから、

(2.23)#\[ \nabla \hat{L}_{\mathcal{D}}(\pmb{w}) = -2 \pmb{X}^\top \pmb{\epsilon} = 2 \pmb{X}^\top (\hat{\pmb{y}} - \pmb{y}) = 2 \pmb{X}^\top (\pmb{X}\pmb{w} - \pmb{y}) \]

目的関数を最小化するパラメータベクトル\(\pmb{w}\)を求めるため、これを\(\pmb{0}\)とおき、\(\pmb{w}\)について解くと、

(2.24)#\[\begin{align} 2\pmb{X}^\top (\pmb{X}\pmb{w} - \pmb{y}) &= \pmb{0} \\ \pmb{X}^\top \pmb{X}\pmb{w} - \pmb{X}^\top \pmb{y} &= \pmb{0} \\ \pmb{X}^\top \pmb{X}\pmb{w} &= \pmb{X}^\top \pmb{y} \\ \pmb{w} &= (\pmb{X}^\top \pmb{X})^{-1}\pmb{X}^\top \pmb{y} \end{align}\]

したがって、学習データに対応する計画行列\(\pmb{X}\)と目的変数ベクトル\(\pmb{y}\)が与えられると、平均二乗残差を最小にするパラメータベクトル\(\pmb{w}\)は行列演算で解析的に(閉じた式で)求めることができる。なお、\((\pmb{X}^\top \pmb{X})^{-1}\pmb{X}^\top\)を求める部分は、ムーア・ペンローズの擬似逆行列(一般化逆行列)でも十分であるが、ここでは詳しく説明しない。

まとめとして、\(\pmb{w}\)を求める式を再掲する。

重回帰のパラメータを求める式

(2.25)#\[ \begin{align} \pmb{w} &= (\pmb{X}^\top \pmb{X})^{-1}\pmb{X}^\top \pmb{y} \end{align} \]

なお、導出の途中で出てきた以下の等式は正規方程式(normal equation)と呼ばれる。

正規方程式

(2.26)#\[ \pmb{X}^\top \pmb{X}\pmb{w} = \pmb{X}^\top \pmb{y} \]

2.5. 目的関数の微分(別の導出)#

(2.19)の目的関数をベクトルで記述したまま、\(\pmb{w}\)で微分する導出も紹介しておく。まず、目的関数を展開する。

(2.27)#\[\begin{align} \hat{L}_{\mathcal{D}}(\pmb{w}) &= \lVert \pmb{y} - \pmb{X}\pmb{w} \rVert^2 \\ &= (\pmb{y} - \pmb{X}\pmb{w})^\top(\pmb{y} - \pmb{X}\pmb{w}) \\ &= (\pmb{y}^\top - \pmb{w}^\top \pmb{X}^\top)(\pmb{y} - \pmb{X}\pmb{w}) \\ &= \pmb{y}^\top \pmb{y} - \pmb{y}^\top \pmb{X}\pmb{w} - \pmb{w}^\top \pmb{X}^\top \pmb{y} + \pmb{w}^\top \pmb{X}^\top \pmb{X}\pmb{w} \\ &= \pmb{y}^\top \pmb{y} - (\pmb{X}^\top\pmb{y})^\top \pmb{w} - \pmb{w}^\top (\pmb{X}^\top \pmb{y}) + \pmb{w}^\top \pmb{X}^\top \pmb{X}\pmb{w} && (\because (\pmb{X}^\top\pmb{y})^\top = \pmb{y}^\top (\pmb{X}^\top)^\top)\\ &= \pmb{y}^\top \pmb{y} - 2(\pmb{X}^\top\pmb{y})^\top \pmb{w} + \pmb{w}^\top \pmb{X}^\top \pmb{X}\pmb{w} && (\because \pmb{u}^\top \pmb{v} = \pmb{v}^\top \pmb{u}) \end{align}\]

続いて、目的関数\(\hat{L}_{\mathcal{D}}(\pmb{w})\)を最小化する\(\pmb{w}\)を求めるために、\(\hat{L}_{\mathcal{D}}(\pmb{w})\)\(\pmb{w}\)で偏微分する。

(2.28)#\[\begin{align} \nabla \hat{L}_{\mathcal{D}}(\pmb{w}) &= \frac{\partial}{\partial \pmb{w}} \left(\pmb{y}^\top \pmb{y} - 2(\pmb{X}^\top\pmb{y})^\top \pmb{w} + \pmb{w}^\top \pmb{X}^\top \pmb{X}\pmb{w}\right) \\ &= 0 - \frac{\partial}{\partial \pmb{w}}\left(2(\pmb{X}^\top\pmb{y})^\top \pmb{w}\right) + \frac{\partial}{\partial \pmb{w}}\left(\pmb{w}^\top \pmb{X}^\top \pmb{X}\pmb{w}\right) \end{align}\]

上式の第2項は、

(2.29)#\[\begin{align} \frac{\partial}{\partial \pmb{w}}\left(2(\pmb{X}^\top\pmb{y})^\top \pmb{w}\right) = 2\pmb{X}^\top \pmb{y} \end{align}\]

第3項は、\(\pmb{X}^\top \pmb{X}\)が正方行列(\(d \times d\)行列)であり、\((\pmb{X}^\top \pmb{X})^\top=\pmb{X}^\top (\pmb{X}^\top)^\top = \pmb{X}^\top \pmb{X}\)、すなわち\(\pmb{X}^\top \pmb{X}\)は対称行列であることに注意すると、

(2.30)#\[\begin{align} \frac{\partial}{\partial \pmb{w}}\left(\pmb{w}^\top \pmb{X}^\top \pmb{X}\pmb{w}\right) = \frac{\partial}{\partial \pmb{w}}\left(\pmb{w}^\top (\pmb{X}^\top \pmb{X}) \pmb{w}\right) = \left((\pmb{X}^\top \pmb{X}) + (\pmb{X}^\top \pmb{X})^\top \right) \pmb{w} = 2\pmb{X}^\top \pmb{X}\pmb{w} \end{align}\]

ゆえに、

(2.31)#\[\begin{align} \nabla \hat{L}_{\mathcal{D}}(\pmb{w}) &= - \frac{\partial}{\partial \pmb{w}}\left(2(\pmb{X}^\top\pmb{y})^\top \pmb{w}\right) + \frac{\partial}{\partial \pmb{w}}\left(\pmb{w}^\top \pmb{X}^\top \pmb{X}\pmb{w}\right) \\ &= - 2\pmb{X}^\top \pmb{y} + 2\pmb{X}^\top \pmb{X}\pmb{w} \\ &= 2\pmb{X}^\top (\pmb{X}\pmb{w} - \pmb{y}) \end{align}\]

したがって、式(2.23)で求めた勾配と一致する。

2.6. 目的関数は凸関数#

ここで、重回帰の目的関数は凸関数であることを示しておく。証明の方針は、目的関数の二階微分(ヘッセ行列)が半正定値行列であることを示すことである。

(2.23)より、目的関数の一階微分は、

(2.32)#\[\begin{align} \nabla \hat{L}_{\mathcal{D}}(\pmb{w}) = \begin{pmatrix} \frac{\partial}{\partial w_1} \\ \frac{\partial}{\partial w_2} \\ \vdots \\ \frac{\partial}{\partial w_d} \\ \end{pmatrix} \hat{L}_{\mathcal{D}}(\pmb{w}) = \frac{\partial \hat{L}_{\mathcal{D}}(\pmb{w})}{\partial \pmb{w}} = 2 \pmb{X}^\top (\pmb{X}\pmb{w} - \pmb{y}) = 2 \pmb{X}^\top\pmb{X}\pmb{w} - 2 \pmb{X}^\top\pmb{y} \end{align}\]

さらに、目的関数の二階微分(ヘッセ行列)は、

(2.33)#\[\begin{align} \nabla^2 \hat{L}_{\mathcal{D}}(\pmb{w}) &= \begin{pmatrix} \frac{\partial^2}{\partial^2 w_1} & \frac{\partial^2}{\partial w_1 \partial w_2} & \dots & \frac{\partial^2}{\partial w_1 \partial w_d} \\ \frac{\partial^2}{\partial w_2 \partial w_1} & \frac{\partial^2}{\partial^2 w_2} & \dots & \frac{\partial^2}{\partial w_2 \partial w_d} \\ \vdots & \vdots & \ddots & \vdots \\ \frac{\partial^2}{\partial w_d \partial w_1} & \frac{\partial^2}{\partial w_d \partial w_2} & \dots & \frac{\partial^2}{\partial^2 w_d} \end{pmatrix} \hat{L}_{\mathcal{D}}(\pmb{w}) \\ &= \frac{\partial}{\partial \pmb{w}} (2 \pmb{X}^\top \pmb{X}\pmb{w} - 2 \pmb{X}^\top \pmb{y}) \\ &= 2 (\pmb{X}^\top \pmb{X})^\top \\ &= 2 \pmb{X}^\top \pmb{X} && (\because \pmb{X}^\top \pmb{X} = \pmb{X}^\top (\pmb{X}^\top)^\top) \end{align}\]

ここで、任意の非ゼロベクトル\(\pmb{u}\)に対して、

(2.34)#\[\begin{align} \pmb{u}^\top (2 \pmb{X}^\top \pmb{X}) \pmb{u} = 2 \pmb{u}^\top \pmb{X}^\top \pmb{X} \pmb{u} = 2(\pmb{X}\pmb{u})^\top (\pmb{X} \pmb{u}) = 2 \|\pmb{X} \pmb{u}\|^2 \geq 0 \end{align}\]

であるから、目的関数の二階微分\(2 \pmb{X}^\top \pmb{X}\)は半正定値行列である。ゆえに、重回帰の目的関数は凸関数である。

2.7. 重回帰の性質#

重回帰に関する様々な性質を導出してみよう [Beck, 2009]。重回帰で推定されたパラメータ\(\pmb{w}\)は、式(2.23)\(\pmb{0}\)とおいたものであるので、

(2.35)#\[\begin{split} \begin{align} \nabla \hat{L}_{\mathcal{D}}(\pmb{w}) = -2 \pmb{X}^\top \pmb{\epsilon} &= \pmb{0} \\ \pmb{X}^\top \pmb{\epsilon} &= \pmb{0} \end{align} \end{split}\]

この行列・ベクトルの成分を明示すると、

(2.36)#\[\begin{align} \begin{pmatrix} X_{1,1} & X_{2,1} & \dots & X_{N,1} \\ X_{1,2} & X_{2,2} & \dots & X_{N,2} \\ \vdots & \vdots & \ddots & \vdots \\ X_{1,d} & X_{2,d} & \dots & X_{N,d} \end{pmatrix} \begin{pmatrix} \epsilon_1 \\ \epsilon_2 \\ \vdots \\ \epsilon_N \end{pmatrix} = \begin{pmatrix} \sum_{i=1}^N X_{i,1} \epsilon_i \\ \sum_{i=1}^N X_{i,2} \epsilon_i \\ \vdots \\ \sum_{i=1}^N X_{i,d} \epsilon_i \end{pmatrix} = \begin{pmatrix} 0 \\ 0 \\ \vdots \\ 0 \end{pmatrix} \end{align}\]

ゆえに、任意の整数\(j \in \{1, 2, \dots, d\}\)に関して、

(2.37)#\[ \sum_{i=1}^N X_{i,j} \epsilon_i = 0 \]

これは、\(\pmb{w}\)を求めるときに、式(2.20)の2行目を全ての\(j\)に対して\(0\)としたことからも、明らかである。

ここで、バイアス項(切片)に対応するパラメータを持たせるため、\(\pmb{X}\)の先頭列の要素はすべて\(1\)と仮定する。

(2.38)#\[ \begin{array}{cc} X_{i,1} = 1 & (i \in \{1, 2, \dots, N\}) \end{array} \]

念のため、\(\pmb{X}\)\(\pmb{X}^\top\)を明示的に書くと、

(2.39)#\[\begin{align} \pmb{X} = \begin{pmatrix} 1 & X_{1,2} & \dots & X_{1,d} \\ 1 & X_{2,2} & \dots & X_{2,d} \\ \vdots & \vdots & \ddots & \vdots \\ 1 & X_{1,N} & \dots & X_{N,d} \end{pmatrix}, \pmb{X}^\top = \begin{pmatrix} 1 & 1 & \dots & 1 \\ X_{1,2} & X_{2,2} & \dots & X_{N,2} \\ \vdots & \vdots & \ddots & \vdots \\ X_{1,d} & X_{2,d} & \dots & X_{N,d} \end{pmatrix} \end{align}\]

これにより、\(w_1\)がバイアス項に対応するパラメータとなる。ここから、重回帰の様々な性質を導出してみよう。なお、バイアス項に対応させる列は先頭である必要はなく、その他の場所でも以降の議論の一般性は失われない。

2.7.1. 残差の和および平均は\(0\)#

\(j=1\)として式(2.37)に式(2.38)を代入すると、

(2.40)#\[ \sum_{i=1}^N X_{i,1} \epsilon_i = \sum_{i=1}^N \epsilon_i = 0 \]

ゆえに、残差の和および平均は\(0\)である。

2.7.2. 目的変数の推定値\(\hat{y}\)の平均値は観測値\(y\)の平均に等しい#

残差の定義\(\pmb{\epsilon} = \pmb{y} - \hat{\pmb{y}}\)を明示的に書くと、

(2.41)#\[\begin{align} \begin{pmatrix} \epsilon_1 \\ \epsilon_2 \\ \vdots \\ \epsilon_N \end{pmatrix} = \begin{pmatrix} y_1 \\ y_2 \\ \vdots \\ y_N \end{pmatrix} - \begin{pmatrix} \hat{y}_1 \\ \hat{y}_2 \\ \vdots \\ \hat{y}_N \end{pmatrix} \end{align}\]

両辺の各要素の和をとると、

(2.42)#\[\begin{align} \sum_{i=1}^N \epsilon_i = \sum_{i=1}^N y_i - \sum_{i=1}^N \hat{y}_i \end{align}\]

ここで、式(2.40)より左辺は\(0\)であるから、

(2.43)#\[\begin{align} \sum_{i=1}^N y_i &= \sum_{i=1}^N \hat{y}_i \\ \frac{1}{N}\sum_{i=1}^N y_i &= \frac{1}{N}\sum_{i=1}^N \hat{y}_i \\ \overline{y} &= \overline{\hat{y}} \end{align}\]

ゆえに、目的変数の推定値\(\hat{y}\)の平均値は観測値\(y\)の平均に等しい。

2.7.3. 観測データの重心は回帰平面上に存在する#

説明変数の平均ベクトル\(\overline{\pmb{x}}\)

(2.44)#\[\begin{align} \overline{\pmb{x}} = \frac{1}{N} \sum_{i=1}^N \pmb{x}_i \end{align}\]

に対して、目的変数の推定値を計算すると、

(2.45)#\[\begin{align} \overline{\pmb{x}}^\top \pmb{w} &= \frac{1}{N} \sum_{i=1}^N \pmb{x}_i^\top \pmb{w} \\ &= \frac{1}{N} \sum_{i=1}^N \hat{y}_i \\ &= \overline{\hat{y}} = \overline{y} \end{align}\]

目的変数の平均\(\overline{y}\)が得られる。ゆえに、観測データの重心は回帰平面上に存在する。

2.7.4. 各説明変数と残差には相関がない#

\(j\)番目の説明変数\(\pmb{X}_{:,j}\)を確率変数と見なし、\(\mathcal{X}_j\)と書くことにする。 任意の説明変数\(\mathcal{X}_j\)と残差\(\mathcal{E}\)の共分散を計算する。残差の平均は\(\overline{\epsilon}=0\)であること、式(2.37)と式(2.40)を利用すると、

(2.46)#\[\begin{align} \mathrm{Cov}[\mathcal{X}_j, \mathcal{E}] &= \frac{1}{N}\sum_{i=1}^N (X_{i,j} - \overline{\pmb{X}_{:,j}}) (\epsilon_i - \overline{\epsilon}) \\ &= \frac{1}{N}\sum_{i=1}^N (X_{i,j} - \overline{\pmb{X}_{:,j}}) \epsilon_i \\ &= \frac{1}{N}\sum_{i=1}^N X_{i,j} \epsilon_i - \frac{1}{N}\sum_{i=1}^N \overline{\pmb{X}_{:,j}} \epsilon_i \\ &= \frac{1}{N}\underbrace{\sum_{i=1}^N X_{i,j} \epsilon_i}_{0} - \frac{\overline{\pmb{X}_{:,j}}}{N}\underbrace{\sum_{i=1}^N \epsilon_i}_{0} \\ &= 0 \end{align}\]

2.7.5. 目的変数の推定値と残差には相関がない#

目的変数の推定値\(\hat{\mathcal{Y}}\)と残差\(\mathcal{E}\)の共分散を計算する。残差の平均は\(\overline{\epsilon}=0\)であること、式(2.35)と式(2.40)を利用すると、

(2.47)#\[\begin{split} \begin{align} \mathrm{Cov}[\hat{\mathcal{Y}},\mathcal{E}] &= \frac{1}{N}\sum_{i=1}^N (\hat{y}_i - \overline{\hat{y}}) (\epsilon_i - \overline{\epsilon}) \\ &= \frac{1}{N}\sum_{i=1}^N (\hat{y}_i - \overline{\hat{y}}) \epsilon_i \\ &= \frac{1}{N}\sum_{i=1}^N \hat{y}_i \epsilon_i - \frac{1}{N}\sum_{i=1}^N \overline{\hat{y}} \epsilon_i \\ &= \frac{1}{N}\sum_{i=1}^N \hat{y}_i \epsilon_i - \frac{\overline{\hat{y}}}{N}\underbrace{\sum_{i=1}^N \epsilon_i}_{0} \\ &= \frac{1}{N}\sum_{i=1}^N \hat{y}_i \epsilon_i \\ &= \frac{1}{N}\pmb{\hat{y}}^\top \pmb{\epsilon} \\ &= \frac{1}{N}(\pmb{X} \pmb{w})^\top \pmb{\epsilon} \\ &= \frac{1}{N}\pmb{w}^\top \underbrace{\pmb{X}^\top \pmb{\epsilon}}_{0} \\ &= 0 \end{align} \end{split}\]

であるから、\(\mathrm{Cov}[\hat{\mathcal{Y}},\mathcal{E}] = 0\)が言える。すなわち、目的変数の推定値と残差の共分散が0になることから、目的変数の推定値と残差は無相関(無関係)であることが分かる。

2.8. 決定係数#

目的変数の観測値の分散\(\mathrm{Var}[\mathcal{Y}]\)は目的変数の推定値の分散\(\mathrm{Var}[\hat{\mathcal{Y}}]\)と残差の分散\(\mathrm{Var}[\mathcal{E}]\)の和で表されることを示す。

(2.48)#\[\begin{align} \mathrm{Var}[\mathcal{Y}] &= \frac{1}{N} \sum_{i=1}^N (y_i - \overline{y})^2 \\ &= \frac{1}{N} \sum_{i=1}^N \{(\hat{y}_i + \epsilon_i) - \overline{y}\}^2 \\ &= \frac{1}{N} \sum_{i=1}^N \{(\hat{y}_i - \overline{y}) + \epsilon_i\}^2 \\ &= \frac{1}{N} \sum_{i=1}^N (\hat{y}_i - \overline{y})^2 + \frac{1}{N} \sum_{i=1}^N 2 (\hat{y}_i - \overline{y}) \epsilon_i + \frac{1}{N} \sum_{i=1}^N \epsilon_i^2 \\ &= \mathrm{Var}[\hat{\mathcal{Y}}] + 2 \cdot \frac{1}{N} \underbrace{\sum_{i=1}^N \hat{y}_i \epsilon_i}_{0} - 2\overline{y} \cdot \frac{1}{N} \underbrace{\sum_{i=1}^N \epsilon_i}_{0} + \mathrm{Var}[\mathcal{E}] \\ &= \mathrm{Var}[\hat{\mathcal{Y}}] + \mathrm{Var}[\mathcal{E}] \end{align}\]

なお、以上の式変形における最後から2行目において、第2項は式(2.47)の導出結果、第3項は式(2.40)より\(0\)であることを利用した。

したがって、単回帰のときと同様に、

  • 全変動: 目的変数の観測値の分散\(\mathrm{Var}[\mathcal{Y}]\)(データの変動)

  • 回帰変動: 目的変数の推定値の分散\(\mathrm{Var}[\hat{\mathcal{Y}}]\)(モデルの変動)

  • 残差変動: 残差の分散\(\mathrm{Var}[\mathcal{E}]\)(モデルが説明しきれなかった変動)

すると、先ほどの式は「全変動は回帰変動と残差変動の和で表される」という関係を表している。

データの全変動のうち、回帰によって表現できた回帰変動の割合を決定係数とよび、次式で定義する。

決定係数

(2.49)#\[\begin{align} R^2 = \frac{\mathrm{Var}[\hat{\mathcal{Y}}]}{\mathrm{Var}[\mathcal{Y}]} = \frac{\mathrm{Var}[\mathcal{Y}] - \mathrm{Var}[\mathcal{E}]}{\mathrm{Var}[\mathcal{Y}]} = 1 - \frac{\mathrm{Var}[\mathcal{E}]}{\mathrm{Var}[\mathcal{Y}]} \end{align}\]

決定係数は\([0,1]\)の範囲の値ととり、データがモデル(単回帰の場合は直線)によく当てはまるときは\(1\)に近くなる。

2.9. 重回帰の実施例#

冒頭の最高気温とアイスクリーム・シャーベットの支出額のデータに対して、単回帰を行う例を示す。まず、説明変数(最高気温)をX、目的変数(支出額)をYに格納する。

def plot_graph(X, Y, x, y):
    fig, ax = plt.subplots(dpi=100)
    ax.set_title('最高気温とアイスクリーム・シャーベットの支出額')
    ax.set_xlabel('最高気温の月平均(℃)')
    ax.set_ylabel('支出額(円)')
    ax.set_xlim(0, 35)
    ax.set_ylim(-250, 2000)
    ax.grid()
    ax.scatter(X, Y, marker='.')
    ax.plot(x, y, 'r')

X = np.array([
     9.1, 11.2, 12.3, 18.9, 22.2, 26. , 30.9, 31.2, 28.8, 23. , 18.3,
    11.1,  8.3,  9.1, 12.5, 18.5, 23.6, 24.8, 30.1, 33.1, 29.8, 23. ,
    16.3, 11.2,  9.6, 10.3, 16.4, 19.2, 24.1, 26.5, 31.4, 33.2, 28.8,
    23. , 17.4, 12.1, 10.6,  9.8, 14.5, 19.6, 24.7, 26.9, 30.5, 31.2,
    26.9, 23. , 17.4, 11. , 10.4, 10.4, 15.5, 19.3, 26.4, 26.4, 30.1,
    30.5, 26.4, 22.7, 17.8, 13.4, 10.6, 12.2, 14.9, 20.3, 25.2, 26.3,
    29.7, 31.6, 27.7, 22.6, 15.5, 13.8, 10.8, 12.1, 13.4, 19.9, 25.1,
    26.4, 31.8, 30.4, 26.8, 20.1, 16.6, 11.1,  9.4, 10.1, 16.9, 22.1,
    24.6, 26.6, 32.7, 32.5, 26.6, 23. , 17.7, 12.1, 10.3, 11.6, 15.4,
    19. , 25.3, 25.8, 27.5, 32.8, 29.4, 23.3, 17.7, 12.6, 11.1, 13.3,
    16. , 18.2, 24. , 27.5, 27.7, 34.1, 28.1, 21.4, 18.6, 12.3])

Y = np.array([
    463.,  360.,  380.,  584.,  763.,  886., 1168., 1325.,  847.,
    542.,  441.,  499.,  363.,  327.,  414.,  545.,  726.,  847.,
   1122., 1355.,  916.,  571.,  377.,  465.,  377.,  362.,  518.,
    683.,  838., 1012., 1267., 1464., 1000.,  629.,  448.,  466.,
    404.,  343.,  493.,  575.,  921., 1019., 1149., 1303.,  805.,
    739.,  587.,  561.,  486.,  470.,  564.,  609.,  899.,  946.,
   1295., 1325.,  760.,  667.,  564.,  633.,  478.,  450.,  567.,
    611.,  947.,  962., 1309., 1307.,  930.,  668.,  496.,  650.,
    506.,  423.,  531.,  672.,  871.,  986., 1368., 1319.,  924.,
    716.,  651.,  708.,  609.,  535.,  717.,  890., 1054., 1077.,
   1425., 1378.,  900.,  725.,  554.,  542.,  561.,  459.,  604.,
    745., 1105.,  973., 1263., 1533., 1044.,  821.,  621.,  601.,
    549.,  572.,  711.,  819., 1141., 1350., 1285., 1643., 1133.,
    784.,  682.,  587.])

2.9.1. numpy.polynomialを用いる例#

numpy.polynomialを用いると、線形回帰モデルのパラメータを求めることができる。多項式近似の次元数を3番目の引数に指定する。以下は2次関数にフィッティングする例である。

from numpy.polynomial import Polynomial

W = Polynomial.fit(X, Y, 2)
W.convert()
\[x \mapsto \text{635.7482577359643} - \text{33.143024247238735}\,x + \text{1.7316079619689542}\,x^{2}\]

多項式近似の係数を\(0\)次、\(1\)次、…の順で並べて表示する。

W.convert().coef
array([635.74825774, -33.14302425,   1.73160796])

\(0\)から\(35\)までに等間隔に配置した\(100\)個の数値xに対して、目的変数の推定値y_hatを求め、回帰曲線を描画する。

x = np.linspace(0, 35, 100)
y_hat = W(x)
plot_graph(X, Y, x, y_hat)
../_images/02mra_50_0.png

続いて、3次関数にフィッティングしてみる。

W = Polynomial.fit(X, Y, 3)
W.convert().coef
array([ 2.64628814e+02,  2.92887138e+01, -1.44685888e+00,  5.01164586e-02])
x = np.linspace(0, 35, 100)
y_hat = W(x)
plot_graph(X, Y, x, y_hat)
../_images/02mra_53_0.png

さらに、4次関数にフィッティングした結果を示す。

W = Polynomial.fit(X, Y, 4)
W.convert().coef
array([-3.47044444e+02,  1.67566320e+02, -1.23937689e+01,  4.12591380e-01,
       -4.27358629e-03])
x = np.linspace(0, 35, 100)
y_hat = W(x)
plot_graph(X, Y, x, y_hat)
../_images/02mra_56_0.png

2.9.2. sklearn.linear_model.LinearRegressionを用いる例#

from sklearn.linear_model import LinearRegression
from sklearn.preprocessing import PolynomialFeatures
from sklearn.pipeline import make_pipeline

sklearn.linear_model.LinearRegressionクラスを用いても、線形回帰を容易に実行できる。訓練データにフィッティングさせるには、fit関数を呼び出せばよいが、説明変数のデータとして、事例数を行数、目的変数の数を列数とした行列形式で与えることになっているので、Xの形状を変更する。以下に、3次関数にフィッティングする例を示す。

reg = make_pipeline(PolynomialFeatures(3), LinearRegression())
reg.fit(X.reshape(-1, 1), Y)
Pipeline(steps=[('polynomialfeatures', PolynomialFeatures(degree=3)),
                ('linearregression', LinearRegression())])
In a Jupyter environment, please rerun this cell to show the HTML representation or trust the notebook.
On GitHub, the HTML representation is unable to render, please try loading this page with nbviewer.org.

回帰直線の係数(傾き)を表示する。

reg[1].coef_
array([ 0.        , 29.28871381, -1.44685888,  0.05011646])

回帰直線の切片を表示する。

reg[1].intercept_
264.628814101544

決定係数を求めるため、score関数を呼び出す。

reg.score(X.reshape(-1, 1), Y)
0.8784609747360882

\(0\)から\(35\)までに等間隔に配置した\(100\)個の数値xに対して、目的変数の推定値y_hatを求めるため、predict関数を呼び出し、回帰曲線を描画する。

x = np.linspace(0, 35, 100)
y_hat = reg.predict(x.reshape(-1, 1))
plot_graph(X, Y, x, y_hat)
../_images/02mra_68_0.png

2.10. 確認問題#

データ\(\mathcal{D}_s\)に対して多項式のフィッティングを行いたい。以下の処理を行うプログラムを作成せよ。

D = np.array([[1, 3], [3, 6], [6, 5], [8, 7]])

なお、NumPy, SciPy, scikit-learn, statsmodel等のライブラリには回帰分析を行う便利な関数として以下のようなものがあるが、ここでは使わずに講義中で説明した式をプログラムとして表現すること。

(1) 行列による1次関数のパラメータ推定

\(y = w_0 + w_1 x\)とおき、学習データ\(\mathcal{D}_s\)上の平均二乗残差を最小にする\(\pmb{w} = \begin{pmatrix} w_0 \\ w_1 \end{pmatrix}\)を行列演算により求めよ。

(2) 2次関数による重回帰

\(y = w_0 + w_1 x + w_2 x^2\)とおき、重回帰により平均二乗残差を最小にする\(\pmb{w} = \begin{pmatrix} w_0 \\ w_1 \\ w_2 \end{pmatrix}\)を求めよ。

(3) 回帰曲線の描画

回帰で求めた2次関数をデータ点とともにグラフに描け。

(4) 決定係数

回帰で得られた2次関数の決定係数(\(R^2\))を求めよ。

(5) 3次関数による重回帰

\(y = w_0 + w_1 x + w_2 x^2 + w_3 x^3\)とおき、重回帰により平均二乗残差を最小にする\(\pmb{w} = \begin{pmatrix} w_0 \\ w_1 \\ w_2 \\ w_3 \end{pmatrix}\)を求めよ。

(6) 回帰曲線の描画

回帰で求めた3次関数をデータ点とともにグラフに描け

(7) 決定係数

回帰で求めた3次関数の決定係数(\(R^2\))を求めよ