重回帰

2. 重回帰#

import numpy as np
import matplotlib.pyplot as plt
import japanize_matplotlib

2.1. 重回帰#

これまで、以下のデータ $D_{s}$ を例として単回帰を説明した。

(2.1)#

\begin{array}{r} D_{s} = {(x_{i}, y_{i})}_{i = 1}^{4} = {(1, 3), (3, 6), (6, 5), (8, 7)} \end{array}

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

../_images/9f1410dfaeddd1f134869c812db5467ee07d67d83cca8c5e8cb0bb71125c85e6.png

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

(2.2)#

\begin{array}{r} \hat{y} = w_{0} + w_{1} x + w_{2} x^{2} + \dots + w_{d} x^{d} \end{array}

にフィッティングすることを考えたい。ここで、 $w_{0}, w_{1}, w_{2}, \dots, w_{d} \in 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}$ の線形結合で表されているため、線形モデルと見なすことができる。ただし、説明変数が２個以上であるので、単回帰を重回帰に拡張する必要がある。

一般的に、 $d$ 個の説明変数 $x_{1}, x_{2}, \dots, x_{d}$ に対して、線形モデル、すなわち $d$ 個の重み $w_{1}, w_{2}, \dots, w_{d}$ の線形結合で目的変数の予測値 $\hat{y}$ を計算することを線形回帰と呼ぶ。

(2.3)#

\begin{array}{r} \hat{y} = w_{1} x_{1} + w_{2} x_{2} + \dots + w_{d} x_{d} \end{array}

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

../_images/9b2327d00574643a1c7bdd115761feeca9020d68078d6ea65037cab92ad9be74.png

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

non-linear

2.2. 表記法の定義#

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

(2.4)#

\begin{array}{r} x = (\begin{array}{c} x_{1} \\ x_{2} \\ ⋮ \\ x_{d} \end{array}) \end{array}

また、 $d$ 個のパラメータ（重み） $w_{1}, w_{2}, \dots, w_{d}$ をベクトル $w \in R^{d}$ でまとめて表現する。

(2.5)#

\begin{array}{r} w = (\begin{array}{c} w_{1} \\ w_{2} \\ ⋮ \\ w_{d} \end{array}) \end{array}

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

(2.6)#

\begin{array}{r} \hat{y} = x^{⊤} w = w_{1} x_{1} + w_{2} x_{2} + \dots + w_{d} x_{d} \end{array}

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

(2.7)#

\begin{array}{r} D = {(x_{1}, y_{1}), (x_{2}, y_{2}), \dots, (x_{N}, y_{N})} = {(x_{i}, y_{i})}_{i = 1}^{N} \end{array}

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

(2.8)#

\begin{array}{r} X = (\begin{array}{c} x_{1}^{⊤} \\ x_{2}^{⊤} \\ ⋮ \\ x_{N}^{⊤} \end{array}) = (\begin{array}{c} X_{1, 1} & X_{1, 2} & \dots & X_{1, d} \\ X_{2, 1} & X_{2, 2} & \dots & X_{2, d} \\ ⋮ & ⋮ & ⋱ & ⋮ \\ X_{N, 1} & X_{N, 2} & \dots & X_{N, d} \end{array}) \end{array}

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

また、目的変数 $y_{i}$ を縦に $N$ 個並べたベクトル $y \in R^{N}$ を定義する。

(2.9)#

\begin{array}{r} y = (\begin{array}{c} y_{1} \\ y_{2} \\ ⋮ \\ y_{N} \end{array}) \end{array}

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

(2.10)#

\begin{array}{r} \hat{y} = (\begin{array}{c} {\hat{y}}_{1} \\ {\hat{y}}_{2} \\ ⋮ \\ {\hat{y}}_{N} \end{array}) = (\begin{array}{c} x_{1}^{⊤} w \\ x_{2}^{⊤} w \\ ⋮ \\ x_{N}^{⊤} w \end{array}) = X w \end{array}

2.2.1. 単回帰を表現する#

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

(2.11)#

\begin{array}{r} \begin{matrix} x = (\begin{array}{c} 1 \\ x \end{array}), w = (\begin{array}{c} b \\ a \end{array}) \end{matrix} \end{array}

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

(2.12)#

\begin{array}{r} \hat{y} = x^{⊤} w = b + x a = a x + b \end{array}

したがって、この章で説明する性質（例えば目的関数が凸関数であること等）は単回帰にも当てはまる。なお、冒頭のグラフの例（ $N = 4$ ）の学習データは次のように表現される。

\begin{array}{r} D_{s} = {(x_{i}, y_{i})}_{i = 1}^{4} = {({(\begin{array}{c} 1 & 1 \end{array})}^{⊤}, 3), ({(\begin{array}{c} 1 & 3 \end{array})}^{⊤}, 6), ({(\begin{array}{c} 1 & 6 \end{array})}^{⊤}, 5), ({(\begin{array}{c} 1 & 8 \end{array})}^{⊤}, 7)} \end{array}

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

(2.13)#

\begin{array}{r} X = (\begin{array}{c} 1 & 1 \\ 1 & 3 \\ 1 & 6 \\ 1 & 8 \end{array}), y = (\begin{array}{c} 3 \\ 6 \\ 5 \\ 7 \end{array}) \end{array}

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

(2.14)#

\begin{array}{r} \hat{y} = X w = (\begin{array}{c} 1 & 1 \\ 1 & 3 \\ 1 & 6 \\ 1 & 8 \end{array}) (\begin{array}{c} b \\ a \end{array}) = (\begin{array}{c} a + b \\ 3 a + b \\ 6 a + b \\ 8 a + b \end{array}) \end{array}

2.3. 目的関数#

単回帰のときと同様に、ある事例の目的変数の真の値 $y_{i}$ と推定値 ${\hat{y}}_{i}$ の差を残差 $ϵ_{i}$ とする。

(2.15)#

\begin{array}{r} ϵ_{i} = y_{i} - {\hat{y}}_{i} = y_{i} - x_{i}^{⊤} w \end{array}

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

(2.16)#

\begin{array}{r} \hat{L} = \frac{1}{N} \sum_{i = 1}^{N} ϵ_{i}^{2} = \frac{1}{N} \sum_{i = 1}^{N} (y_{i} - {\hat{y}}_{i})^{2} = \frac{1}{N} \sum_{i = 1}^{N} (y_{i} - x_{i}^{⊤} w)^{2} \end{array}

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

\begin{array}{r} {\hat{L}}_{D} (w) = \sum_{i = 1}^{N} (y_{i} - x_{i}^{⊤} w)^{2} \end{array}

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

(2.17)#

\begin{array}{r} ϵ = (\begin{array}{c} ϵ_{1} \\ ϵ_{2} \\ ⋮ \\ ϵ_{N} \end{array}) \end{array}

で表すことにすると、

(2.18)#

\begin{array}{r} ϵ = y - \hat{y} = y - X w \end{array}

である。

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

重回帰の目的関数（二乗残差和）

(2.19)#

\begin{array}{r} {\hat{L}}_{D} (w) = ‖ ϵ ‖^{2} = ‖ y - X w ‖^{2} = \sum_{i = 1}^{N} (y_{i} - x_{i}^{⊤} w)^{2} \end{array}

と表現できる。

2.4. 目的関数の最小化#

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

(2.20)#

\begin{array}{r} \begin{aligned} \frac{\partial {\hat{L}}_{D} (w)}{\partial w_{j}} & = \sum_{i = 1}^{N} 2 \cdot \underset{= ϵ_{i}}{\underset{⏟}{(y_{i} - x_{i}^{⊤} w)}} \cdot {- \underset{= X_{i, j}}{\underset{⏟}{(x_{i})_{j}}}} \\ = - 2 \sum_{i = 1}^{N} X_{i, j} ϵ_{i} \\ = - 2 (\begin{array}{c} X_{1, j} & X_{2, j} & \dots & X_{N, j} \end{array}) (\begin{array}{c} ϵ_{1} \\ ϵ_{2} \\ ⋮ \\ ϵ_{N} \end{array}) \\ = - 2 (X^{⊤})_{j} ϵ \\ = - 2 (X^{⊤} ϵ)_{j} \end{aligned} \end{array}

なお、最後から２番目の式変形は、計画行列 $X$ の転置 $X^{⊤}$ 、

(2.21)#

\begin{array}{r} X^{⊤} = (\begin{array}{c} X_{1, 1} & X_{2, 1} & \dots & X_{N, 1} \\ X_{1, 2} & X_{2, 2} & \dots & X_{N, 2} \\ ⋮ & ⋮ & ⋱ & ⋮ \\ X_{1, j} & X_{2, j} & \dots & X_{N, j} \\ ⋮ & ⋮ & ⋱ & ⋮ \\ X_{1, d} & X_{2, d} & \dots & X_{N, d} \end{array}) \end{array}

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

(2.22)#

\begin{array}{r} X^{⊤} ϵ = (\begin{array}{c} (X^{⊤})_{1} ϵ \\ (X^{⊤})_{2} ϵ \\ ⋮ \\ (X^{⊤})_{j} ϵ \\ ⋮ \\ (X^{⊤})_{d} ϵ \end{array}) \end{array}

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

\begin{array}{r} \nabla {\hat{L}}_{D} (w) = (\begin{array}{c} \frac{\partial}{\partial w_{1}} \\ \frac{\partial}{\partial w_{2}} \\ ⋮ \\ \frac{\partial}{\partial w_{d}} \end{array}) {\hat{L}}_{D} (w) = - 2 X^{⊤} ϵ \end{array}

さらに、式(2.18)より $ϵ = y - X w$ であるから、

(2.23)#

\nabla {\hat{L}}_{D} (w) = - 2 X^{⊤} ϵ = 2 X^{⊤} (\hat{y} - y) = 2 X^{⊤} (X w - y)

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

(2.24)#

\begin{aligned} 2 X^{⊤} (X w - y) & = 0 \\ X^{⊤} X w - X^{⊤} y & = 0 \\ X^{⊤} X w & = X^{⊤} y \\ w & = (X^{⊤} X)^{- 1} X^{⊤} y \end{aligned}

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

まとめとして、 $w$ を求める式を再掲する。

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

(2.25)#

\begin{aligned} w & = (X^{⊤} X)^{- 1} X^{⊤} y \end{aligned}

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

正規方程式

(2.26)#

X^{⊤} X w = X^{⊤} y

2.5. 目的関数の微分（別の導出）#

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

(2.27)#

\begin{aligned} {\hat{L}}_{D} (w) & = ‖ y - X w ‖^{2} \\ = (y - X w)^{⊤} (y - X w) \\ = (y^{⊤} - w^{⊤} X^{⊤}) (y - X w) \\ = y^{⊤} y - y^{⊤} X w - w^{⊤} X^{⊤} y + w^{⊤} X^{⊤} X w \\ = y^{⊤} y - (X^{⊤} y)^{⊤} w - w^{⊤} (X^{⊤} y) + w^{⊤} X^{⊤} X w & (∵ (X^{⊤} y)^{⊤} = y^{⊤} (X^{⊤})^{⊤}) \\ = y^{⊤} y - 2 (X^{⊤} y)^{⊤} w + w^{⊤} X^{⊤} X w & (∵ u^{⊤} v = v^{⊤} u) \end{aligned}

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

(2.28)#

\begin{aligned} \nabla {\hat{L}}_{D} (w) & = \frac{\partial}{\partial w} (y^{⊤} y - 2 (X^{⊤} y)^{⊤} w + w^{⊤} X^{⊤} X w) \\ = 0 - \frac{\partial}{\partial w} (2 (X^{⊤} y)^{⊤} w) + \frac{\partial}{\partial w} (w^{⊤} X^{⊤} X w) \end{aligned}

上式の第2項は、

(2.29)#

\begin{array}{r} \frac{\partial}{\partial w} (2 (X^{⊤} y)^{⊤} w) = 2 X^{⊤} y \end{array}

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

(2.30)#

\begin{array}{r} \frac{\partial}{\partial w} (w^{⊤} X^{⊤} X w) = \frac{\partial}{\partial w} (w^{⊤} (X^{⊤} X) w) = ((X^{⊤} X) + (X^{⊤} X)^{⊤}) w = 2 X^{⊤} X w \end{array}

ゆえに、

(2.31)#

\begin{aligned} \nabla {\hat{L}}_{D} (w) & = - \frac{\partial}{\partial w} (2 (X^{⊤} y)^{⊤} w) + \frac{\partial}{\partial w} (w^{⊤} X^{⊤} X w) \\ = - 2 X^{⊤} y + 2 X^{⊤} X w \\ = 2 X^{⊤} (X w - y) \end{aligned}

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

2.6. 目的関数は凸関数#

ここで、重回帰の目的関数は凸関数であることを示しておく。証明の方針は、目的関数の二階微分（ヘッセ行列）が半正定値行列であることを示すことであるが、その方針で証明ができる背景を簡単に説明する。

2.6.1. 凸関数であるための条件#

まずは重回帰の問題設定を忘れ、凸関数に関する一般的な説明を行う。

凸関数の定義

$d$ 次元ベクトル $x \in R^{d}$ を引数にとり、スカラーを返す関数 $f : R^{d} ⟼ R$ がある。任意の $u, v \in R^{d}$ と $0 \leq λ \leq 1$ を満たす任意の $λ$ に対して、

(2.32)#

\begin{array}{r} f (λ u + (1 - λ) v) \leq λ f (u) + (1 - λ) f (v) \end{array}

が成り立つとき、 $f$ を凸関数と呼ぶ。

この定義は「関数 $f$ 上の任意の2点間の曲線は、その2点を結ぶ直線（弦）よりも下にある」ことを表現している（厳密には曲線が弦が同じでも凸関数と呼べる）。 $f$ を1変数関数として、この凸関数の定義を可視化したのが下のグラフである。点 $(u, f (u)), (v, f (v))$ を結ぶ弦は、関数 $f$ の曲線よりも上に位置している。

../_images/93557e8a583ff94144d5628b8661a8f5064f8557e53cc4767db94ff207c97fc9.png

一階微分と凸関数に関する定理

関数 $f : R^{d} ⟼ R$ が一階微分可能であるとする。任意の $x, z \in R^{d}$ に対して、

(2.33)#

\begin{array}{r} f (z) \geq f (x) + \nabla f (x)^{⊤} (z - x) \end{array}

が成り立つとき、 $f$ は凸関数である。

【証明】任意の $u, v \in R^{d}$ と $0 \leq λ \leq 1$ を満たす任意の $λ$ を用いて、 $x = λ u + (1 - λ) v$ と表す。式(2.33)の不等式を $z = u, v$ として適用すると、

\begin{aligned} f (u) \geq & f (x) + \nabla f (x)^{⊤} (u - x) = f (x) + (1 - λ) \nabla f (x)^{⊤} (u - v) \\ f (v) \geq & f (x) + \nabla f (x)^{⊤} (v - x) = f (x) + λ \nabla f (x)^{⊤} (v - u) \end{aligned}

これらの不等式に基づいて $λ f (u) + (1 - λ) f (v)$ を計算すると、

\begin{aligned} λ f (u) + (1 - λ) f (v) & \geq λ {f (x) + (1 - λ) \nabla f (x)^{⊤} (u - v)} + (1 - λ) {f (x) + λ \nabla f (x)^{⊤} (v - u)} \\ = f (x) + λ (1 - λ) \nabla f (x)^{⊤} (u - v) + λ (1 - λ) \nabla f (x)^{⊤} (v - u) \\ = f (x) + λ (1 - λ) \nabla f (x)^{⊤} (u - v + v - u) \\ = f (x) \\ = f (λ u + (1 - λ) v) \end{aligned}

である（最終行の変形では $x$ の $u$ と $v$ による表示を代入した）。この不等式は凸関数の定義（式(2.32)）そのものであり、 $f$ は凸関数であることが示された。なお、ここでは省略するが、この逆を証明して必要十分条件を示すことができる。

この定理は、「関数 $f$ 上の任意の点の接線は $f$ の曲線よりも下にある」ことを表現している（厳密には曲線が接線が同じ点で接することがある）。 $f$ が1変数関数のとき、ある $x, z$ に対してこの定理が成り立つ様子を図示した。

../_images/28d32db733db6056490a239d3a08582a0917ce7a978c1bc3d0fcea71b3ed0db7.png

二階微分と凸関数に関する定理

関数 $f : R^{d} ⟼ R$ が二階微分可能であるとする。 $f$ の二階微分であるヘッセ行列、

\begin{array}{r} \nabla^{2} f (x) = (\begin{array}{c} \frac{\partial^{2} f}{\partial^{2} x_{1}} & \frac{\partial^{2} f}{\partial x_{1} \partial x_{2}} & \dots & \frac{\partial^{2} f}{\partial x_{1} \partial x_{d}} \\ \frac{\partial^{2} f}{\partial x_{2} x_{1}} & \frac{\partial^{2} f}{\partial^{2} x_{2}} & \dots & \frac{\partial^{2} f}{\partial x_{2} \partial x_{d}} \\ ⋮ & ⋮ & ⋱ & ⋮ \\ \frac{\partial^{2} f}{\partial x_{d} x_{1}} & \frac{\partial^{2} f}{\partial x_{d} x_{2}} & \dots & \frac{\partial^{2} f}{\partial^{2} x_{d}} \end{array}) \in R^{d \times d} \end{array}

が半正定値行列であるとき、すなわち、任意のベクトル $u$ に対して、

(2.34)#

\begin{array}{r} u^{⊤} \nabla^{2} f (x) u \geq 0 \end{array}

が成り立つとき、 $f$ は凸関数である。

【証明】テイラーの定理より、任意の $x, z \in R^{d}$ に対して、 $0 < λ < 1$ を満たす $λ$ が存在し、

\begin{array}{r} f (z) = f (x) + \nabla f (x)^{⊤} (z - x) + \frac{1}{2} (z - x)^{⊤} \nabla^{2} f (x + λ (z - x)) (z - x) \end{array}

を満たす。仮定より、 $f$ のヘッセ行列は半正定値行列であることから、上式の第3項は常に非負である（ $u = z - x$ とおけばよい）。ゆえに、

\begin{array}{r} f (z) \geq f (x) + \nabla f (x)^{⊤} (z - x) \end{array}

であり、一階微分と凸関数に関する定理（式(2.33)）から、 $f$ は凸関数であることが示された。なお、ここでは省略するが、この逆を証明して必要十分条件を示すことができる。

$f$ が1変数関数のとき、式(2.34)の条件は $f^{″} (x) \geq 0$ となる。これは、 $f$ の傾きが単調増加する（厳密には減少しない）ことを表している。

2.6.2. 目的関数が凸関数であることの証明#

さて、ヘッセ行列の半正定値性から凸関数であることを示せることが分かったので、変数を $x \to w$ に、目的関数を $f (x) \to {\hat{L}}_{D} (w)$ に置換して考察する。式(2.23)より、目的関数の一階微分は、

(2.35)#

\begin{array}{r} \nabla {\hat{L}}_{D} (w) = (\begin{array}{c} \frac{\partial}{\partial w_{1}} \\ \frac{\partial}{\partial w_{2}} \\ ⋮ \\ \frac{\partial}{\partial w_{d}} \end{array}) {\hat{L}}_{D} (w) = \frac{\partial {\hat{L}}_{D} (w)}{\partial w} = 2 X^{⊤} (X w - y) = 2 X^{⊤} X w - 2 X^{⊤} y \end{array}

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

(2.36)#

\begin{aligned} \nabla^{2} {\hat{L}}_{D} (w) & = (\begin{array}{c} \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}} \\ ⋮ & ⋮ & ⋱ & ⋮ \\ \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{array}) {\hat{L}}_{D} (w) \\ = \frac{\partial}{\partial w} (2 X^{⊤} X w - 2 X^{⊤} y) \\ = 2 (X^{⊤} X)^{⊤} \\ = 2 X^{⊤} X & (∵ (X^{⊤} X) は対称行列) \end{aligned}

ここで、任意のベクトル $u \in R^{d}$ に対して、

(2.37)#

\begin{array}{r} u^{⊤} (2 X^{⊤} X) u = 2 u^{⊤} X^{⊤} X u = 2 (X u)^{⊤} (X u) = 2 ‖ X u ‖^{2} \geq 0 \end{array}

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

なお、凸関数の定義を用いて重回帰の目的関数が凸であることを証明する方法もあるが、ここでは省略する。

2.7. 重回帰の性質#

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

(2.38)#

\begin{array}{r} \begin{aligned} \nabla {\hat{L}}_{D} (w) = - 2 X^{⊤} ϵ & = 0 \\ X^{⊤} ϵ & = 0 \end{aligned} \end{array}

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

(2.39)#

\begin{array}{r} (\begin{array}{c} X_{1, 1} & X_{2, 1} & \dots & X_{N, 1} \\ X_{1, 2} & X_{2, 2} & \dots & X_{N, 2} \\ ⋮ & ⋮ & ⋱ & ⋮ \\ X_{1, d} & X_{2, d} & \dots & X_{N, d} \end{array}) (\begin{array}{c} ϵ_{1} \\ ϵ_{2} \\ ⋮ \\ ϵ_{N} \end{array}) = (\begin{array}{c} \sum_{i = 1}^{N} X_{i, 1} ϵ_{i} \\ \sum_{i = 1}^{N} X_{i, 2} ϵ_{i} \\ ⋮ \\ \sum_{i = 1}^{N} X_{i, d} ϵ_{i} \end{array}) = (\begin{array}{c} 0 \\ 0 \\ ⋮ \\ 0 \end{array}) \end{array}

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

(2.40)#

\sum_{i = 1}^{N} X_{i, j} ϵ_{i} = 0

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

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

(2.41)#

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

念のため、 $X$ と $X^{⊤}$ を明示的に書くと、

(2.42)#

\begin{array}{r} X = (\begin{array}{c} 1 & X_{1, 2} & \dots & X_{1, d} \\ 1 & X_{2, 2} & \dots & X_{2, d} \\ ⋮ & ⋮ & ⋱ & ⋮ \\ 1 & X_{1, N} & \dots & X_{N, d} \end{array}), X^{⊤} = (\begin{array}{c} 1 & 1 & \dots & 1 \\ X_{1, 2} & X_{2, 2} & \dots & X_{N, 2} \\ ⋮ & ⋮ & ⋱ & ⋮ \\ X_{1, d} & X_{2, d} & \dots & X_{N, d} \end{array}) \end{array}

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

2.7.1. 残差の和および平均は $0$ #

$j = 1$ として式(2.40)に式(2.41)を代入すると、

(2.43)#

\sum_{i = 1}^{N} X_{i, 1} ϵ_{i} = \sum_{i = 1}^{N} ϵ_{i} = 0

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

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

残差の定義 $ϵ = y - \hat{y}$ を明示的に書くと、

(2.44)#

\begin{array}{r} (\begin{array}{c} ϵ_{1} \\ ϵ_{2} \\ ⋮ \\ ϵ_{N} \end{array}) = (\begin{array}{c} y_{1} \\ y_{2} \\ ⋮ \\ y_{N} \end{array}) - (\begin{array}{c} {\hat{y}}_{1} \\ {\hat{y}}_{2} \\ ⋮ \\ {\hat{y}}_{N} \end{array}) \end{array}

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

(2.45)#

\begin{array}{r} \sum_{i = 1}^{N} ϵ_{i} = \sum_{i = 1}^{N} y_{i} - \sum_{i = 1}^{N} {\hat{y}}_{i} \end{array}

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

(2.46)#

\begin{aligned} \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} \\ \overset{―}{y} & = \overset{―}{\hat{y}} \end{aligned}

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

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

説明変数の平均ベクトル $\overset{―}{x}$ 、

(2.47)#

\begin{array}{r} \overset{―}{x} = \frac{1}{N} \sum_{i = 1}^{N} x_{i} \end{array}

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

(2.48)#

\begin{aligned} {\overset{―}{x}}^{⊤} w & = \frac{1}{N} \sum_{i = 1}^{N} x_{i}^{⊤} w \\ = \frac{1}{N} \sum_{i = 1}^{N} {\hat{y}}_{i} \\ = \overset{―}{\hat{y}} = \overset{―}{y} \end{aligned}

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

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

$j$ 番目の説明変数 $X_{:, j}$ を確率変数と見なし、 $X_{j}$ と書くことにする。任意の説明変数 $X_{j}$ と残差 $E$ の共分散を計算する。残差の平均は $\overset{―}{ϵ} = 0$ であること、式(2.40)と式(2.43)を利用すると、

(2.49)#

\begin{aligned} Cov [X_{j}, E] & = \frac{1}{N} \sum_{i = 1}^{N} (X_{i, j} - \overset{―}{X_{:, j}}) (ϵ_{i} - \overset{―}{ϵ}) \\ = \frac{1}{N} \sum_{i = 1}^{N} (X_{i, j} - \overset{―}{X_{:, j}}) ϵ_{i} \\ = \frac{1}{N} \sum_{i = 1}^{N} X_{i, j} ϵ_{i} - \frac{1}{N} \sum_{i = 1}^{N} \overset{―}{X_{:, j}} ϵ_{i} \\ = \frac{1}{N} \underset{0}{\underset{⏟}{\sum_{i = 1}^{N} X_{i, j} ϵ_{i}}} - \frac{\overset{―}{X_{:, j}}}{N} \underset{0}{\underset{⏟}{\sum_{i = 1}^{N} ϵ_{i}}} \\ = 0 \end{aligned}

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

目的変数の推定値 $\hat{Y}$ と残差 $E$ の共分散を計算する。残差の平均は $\overset{―}{ϵ} = 0$ であること、式(2.38)と式(2.43)を利用すると、

(2.50)#

\begin{array}{r} \begin{aligned} Cov [\hat{Y}, E] & = \frac{1}{N} \sum_{i = 1}^{N} ({\hat{y}}_{i} - \overset{―}{\hat{y}}) (ϵ_{i} - \overset{―}{ϵ}) \\ = \frac{1}{N} \sum_{i = 1}^{N} ({\hat{y}}_{i} - \overset{―}{\hat{y}}) ϵ_{i} \\ = \frac{1}{N} \sum_{i = 1}^{N} {\hat{y}}_{i} ϵ_{i} - \frac{1}{N} \sum_{i = 1}^{N} \overset{―}{\hat{y}} ϵ_{i} \\ = \frac{1}{N} \sum_{i = 1}^{N} {\hat{y}}_{i} ϵ_{i} - \frac{\overset{―}{\hat{y}}}{N} \underset{0}{\underset{⏟}{\sum_{i = 1}^{N} ϵ_{i}}} \\ = \frac{1}{N} \sum_{i = 1}^{N} {\hat{y}}_{i} ϵ_{i} \\ = \frac{1}{N} {\hat{y}}^{⊤} ϵ \\ = \frac{1}{N} (X w)^{⊤} ϵ \\ = \frac{1}{N} w^{⊤} \underset{0}{\underset{⏟}{X^{⊤} ϵ}} \\ = 0 \end{aligned} \end{array}

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

2.8. 決定係数#

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

(2.51)#

\begin{aligned} Var [Y] & = \frac{1}{N} \sum_{i = 1}^{N} (y_{i} - \overset{―}{y})^{2} \\ = \frac{1}{N} \sum_{i = 1}^{N} {({\hat{y}}_{i} + ϵ_{i}) - \overset{―}{y}}^{2} \\ = \frac{1}{N} \sum_{i = 1}^{N} {({\hat{y}}_{i} - \overset{―}{y}) + ϵ_{i}}^{2} \\ = \frac{1}{N} \sum_{i = 1}^{N} ({\hat{y}}_{i} - \overset{―}{y})^{2} + \frac{1}{N} \sum_{i = 1}^{N} 2 ({\hat{y}}_{i} - \overset{―}{y}) ϵ_{i} + \frac{1}{N} \sum_{i = 1}^{N} ϵ_{i}^{2} \\ = Var [\hat{Y}] + 2 \cdot \frac{1}{N} \underset{0}{\underset{⏟}{\sum_{i = 1}^{N} {\hat{y}}_{i} ϵ_{i}}} - 2 \overset{―}{y} \cdot \frac{1}{N} \underset{0}{\underset{⏟}{\sum_{i = 1}^{N} ϵ_{i}}} + Var [E] \\ = Var [\hat{Y}] + Var [E] \end{aligned}

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

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

全変動: 目的変数の観測値の分散 $Var [Y]$ （データの変動）
回帰変動: 目的変数の推定値の分散 $Var [\hat{Y}]$ （モデルの変動）
残差変動: 残差の分散 $Var [E]$ （モデルが説明しきれなかった変動）

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

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

決定係数

(2.52)#

\begin{array}{r} R^{2} = \frac{Var [\hat{Y}]}{Var [Y]} = \frac{Var [Y] - Var [E]}{Var [Y]} = 1 - \frac{Var [E]}{Var [Y]} \end{array}

決定係数は $[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を用いると、線形回帰モデルのパラメータを求めることができる。多項式近似の次元数を３番目の引数に指定する。以下は2次関数にフィッティングする例である。

from numpy.polynomial import Polynomial

W = Polynomial.fit(X, Y, 2)
W.convert()

x \mapsto 635.7482577359643 - 33.143024247238735 x + 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/29cf1f4af903fac1d996284bc1693bdd3a45e4341bde4da37de3dbc413c1ab0d.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/1b28d94693d72a595d27ef28f7ad2944fcfa45b2405954c0ea050d84da3bca3d.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/e9b97e77154cc1690cbdd5a68310f9e0ef9d7d92548436838002688893a2e4d9.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)

2.10. 確認問題#

データ $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$ とおき、学習データ $D_{s}$ 上の平均二乗残差を最小にする $w = (\begin{matrix} w_{0} \\ w_{1} \end{matrix})$ を行列演算により求めよ。

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

$y = w_{0} + w_{1} x + w_{2} x^{2}$ とおき、重回帰により平均二乗残差を最小にする $w = (\begin{matrix} w_{0} \\ w_{1} \\ w_{2} \end{matrix})$ を求めよ。

(3) 回帰曲線の描画

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

(4) 決定係数

回帰で得られた2次関数の決定係数（ $R^{2}$ ）を求めよ。

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

$y = w_{0} + w_{1} x + w_{2} x^{2} + w_{3} x^{3}$ とおき、重回帰により平均二乗残差を最小にする $w = (\begin{matrix} w_{0} \\ w_{1} \\ w_{2} \\ w_{3} \end{matrix})$ を求めよ。

(6) 回帰曲線の描画

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

(7) 決定係数

回帰で求めた3次関数の決定係数（ $R^{2}$ ）を求めよ

重回帰

目次