線形回帰は，ある目的データを別のいくつかの観測データの線形和として表すことで，両者の関係性を説明し，新たに得られた観測データを元に未知の目的データを推定する用途で用いるモデルです．
線形単回帰モデル
線形単回帰モデルとは，2つのデータ変数 (x x x y y y
統計の基礎 で既にふれているように，その基本は中学校あたりで学習済みです．
Pythonコード
X, y, coef = make_regression( n_samples= 100 , n_features= 1 , noise= 10 , coef= True , random_state= 42 )
model = LinearRegression( )
model. fit( X, y)
y_pred = model. predict( X)
plt. figure( figsize= ( 8 , 6 ) )
plt. scatter( X, y, color= 'blue' , label= 'Data' )
plt. plot( X, y_pred, color= 'red' , label= 'Regression Line' )
plt. title( 'Linear Regression' , fontsize= 20 )
plt. xlabel( 'Feature' , fontsize= 15 )
plt. ylabel( 'Target' , fontsize= 15 )
plt. legend( )
plt. grid( True )
plt. savefig( '../figures/linear.png' )
この図では，てきとうに作成したデータ点に対して単回帰を行い，得られた線 (回帰線) を赤で引いてあります．良い感じですね．たしかにそれっぽいところに線が引かれています．
概要はこの程度で十分なので，実際にモデルの中身を見ていきます．
まず，単回帰の場合は説明変数 x x x 目的変数 y y y
y = β 0 + β 1 x \begin{align}
y = \beta_0 + \beta_1 x
\end{align} y = β 0 + β 1 x
というシンプルな関係性を仮定します．ここで β 0 , 1 \beta_{0,1} β 0 , 1 β 1 \beta_1 β 1 回帰係数 といいます．発想は初等数学でやる y = a x + b y=ax+b y = a x + b ( β 0 , β 1 ) (\beta_0, \beta_1) ( β 0 , β 1 )
ただ，実際には測定の誤差や，今回考えていない他の要因 (x 2 , 3 , . . . x_{2,3,...} x 2 , 3 , ... ϵ \epsilon ϵ
従って線形単回帰モデル はデータの組 (標本) の総数を n n n
y i = β 0 + β 1 x i + ϵ i \begin{align}
y_i = \beta_0 + \beta_1 x_i + \epsilon_i
\end{align} y i = β 0 + β 1 x i + ϵ i
という形になります．ここで, i = 1 , . . . , n i.i.d. i = 1,...,n \quad \text{i.i.d.} i = 1 , ... , n i.i.d.
言葉にすると，得られた i i i y i y_i y i x i x_i x i β 1 \beta_1 β 1 β 0 \beta_0 β 0 ϵ i \epsilon_i ϵ i ，というモデルです．
単回帰係数の推定
では次に，得られた線形単回帰モデルの β 0 , 1 \beta_{0,1} β 0 , 1
最小二乗法
最小二乗法では，
「データに良く当てはまるモデルは，誤差 ϵ i \epsilon_i ϵ i i i i
と考えます．尤もらしいですね．ちなみに，ϵ i \epsilon_i ϵ i
ϵ 1 , . . . , ϵ n i.i.d. E ( ϵ i ) = 0 V ( ϵ i ) = σ 2 \begin{align}
\epsilon_1,...,\epsilon_n \quad \text{i.i.d.}\\
E(\epsilon_i)=0\\
V(\epsilon_i)=\sigma^2
\end{align} ϵ 1 , ... , ϵ n i.i.d. E ( ϵ i ) = 0 V ( ϵ i ) = σ 2
誤差はデータごとに独立で，平均0，分散 σ 2 \sigma^2 σ 2 σ 2 \sigma^2 σ 2
これにより，同時に
y 1 , . . . , y n i.i.d. E ( y i ) = β 0 + β 1 x i V ( y i ) = σ 2 \begin{align}
y_1,...,y_n \quad \text{i.i.d.}\\
E(y_i)=\beta_0 + \beta_1 x_i\\
V(y_i)=\sigma^2
\end{align} y 1 , ... , y n i.i.d. E ( y i ) = β 0 + β 1 x i V ( y i ) = σ 2
となります．
では本題です．
最小二乗法とは，予測値 y i ^ \hat{y_i} y i ^ y i y_i y i を意味し，今回のモデルについて考えると
arg min β 0 , β 1 S = ∑ i = 1 n ϵ i 2 = ∑ i = 1 n { y i − ( β 0 + β 1 x i ) } 2 \argmin_{\beta_0,\beta_1} S = \sum_{i=1}^n \epsilon_i^2 = \sum_{i=1}^n \{ y_i - (\beta_0 + \beta_1 x_i) \}^2 β 0 , β 1 arg min S = i = 1 ∑ n ϵ i 2 = i = 1 ∑ n { y i − ( β 0 + β 1 x i ) } 2
という問題になります．
arg min β 0 , β 1 S \argmin_{\beta_0,\beta_1} S arg min β 0 , β 1 S S S S β 0 , β 1 \beta_0, \beta_1 β 0 , β 1 S S S β 0 , β 1 \beta_0, \beta_1 β 0 , β 1 S S S y i y_i y i x i x_i x i y i ^ \hat{y_i} y i ^ i i i
arg min β 0 , β 1 S \argmin_{\beta_0,\beta_1} S arg min β 0 , β 1 S
変数それぞれについて，関数 S S S β 0 , β 1 \beta_0, \beta_1 β 0 , β 1
∂ S ( β 0 , β 1 ) ∂ β 0 = − 2 ∑ i = 1 n ( y i − β 1 x i − β 0 ) ∂ S ( β 0 , β 1 ) ∂ β 1 = − 2 ∑ i = 1 n x i ( y i − β 1 x i − β 0 ) \begin{align}
\frac{\partial S(\beta_0, \beta_1)}{\partial \beta_0} = -2 \sum_{i=1}^n (y_i - \beta_1 x_i -\beta_0)\\
\frac{\partial S(\beta_0, \beta_1)}{\partial \beta_1} = -2\sum_{i=1}^n x_i(y_i -\beta_1 x_i- \beta_0)
\end{align} ∂ β 0 ∂ S ( β 0 , β 1 ) = − 2 i = 1 ∑ n ( y i − β 1 x i − β 0 ) ∂ β 1 ∂ S ( β 0 , β 1 ) = − 2 i = 1 ∑ n x i ( y i − β 1 x i − β 0 )
となります．で，こいつらの最小値を求めるわけなのでそれぞれを = 0 =0 = 0
β 1 ∑ i = 1 n x i + β 0 n = ∑ i = 1 n y i β 1 ∑ i = 1 n x i 2 + β 0 ∑ i = 1 n x i = ∑ i = 1 n x i y i \begin{align}
\beta_1 \sum_{i=1}^n x_i + \beta_0 n = \sum_{i=1}^n y_i \\
\beta_1 \sum_{i=1}^n x_i^2 + \beta_0 \sum_{i=1}^n x_i = \sum_{i=1}^n x_i y_i
\end{align} β 1 i = 1 ∑ n x i + β 0 n = i = 1 ∑ n y i β 1 i = 1 ∑ n x i 2 + β 0 i = 1 ∑ n x i = i = 1 ∑ n x i y i
が求められます．
1 つ目の式を n n n
β 1 ∑ i = 1 n x i n + β 0 = ∑ i = 1 n y i n β 1 x ˉ + β 0 = y ˉ ∴ β 0 ^ = y ˉ − β 1 ^ x ˉ \frac{\beta_1\sum_{i=1}^n x_i}{n} + \beta_0 = \frac{\sum_{i=1}^n y_i}{n}\\
\beta_1 \bar{x} + \beta_0 = \bar{y}\\
\therefore \hat{\beta_0} = \bar{y} - \hat{\beta_{1}} \bar{x} n β 1 ∑ i = 1 n x i + β 0 = n ∑ i = 1 n y i β 1 x ˉ + β 0 = y ˉ ∴ β 0 ^ = y ˉ − β 1 ^ x ˉ
となります．x ˉ \bar{x} x ˉ x x x n n n
途中式
β 1 ∑ i = 1 n x i 2 + ( − β 1 x ˉ + y ˉ ) ∑ i = 1 n x i = ∑ i = 1 n x i y i β 1 ( ∑ i = 1 n x i 2 n − x ˉ ∑ i = 1 n x i n ) = ∑ i = 1 n x i y i n − y ˉ ∑ i = 1 n x i n a ( x 2 ˉ − x ˉ 2 ) = x y ˉ − y ˉ x ˉ \beta_1\sum_{i=1}^n x_i^2 + (-\beta_1\bar{x} + \bar{y})\sum_{i=1}^n x_i = \sum_{i=1}^n x_i y_i \nonumber\\
\beta_1(\frac{\sum_{i=1}^n x_i^2}{n} - \bar{x}\frac{\sum_{i=1}^n x_i}{n}) = \frac{\sum_{i=1}^n x_i y_i}{n} - \bar{y}\frac{\sum_{i=1}^n x_i}{n}\\
a(\bar{x^2} - \bar{x}^2) = \bar{xy} - \bar{y}\bar{x}\\ β 1 i = 1 ∑ n x i 2 + ( − β 1 x ˉ + y ˉ ) i = 1 ∑ n x i = i = 1 ∑ n x i y i β 1 ( n ∑ i = 1 n x i 2 − x ˉ n ∑ i = 1 n x i ) = n ∑ i = 1 n x i y i − y ˉ n ∑ i = 1 n x i a ( x 2 ˉ − x ˉ 2 ) = x y ˉ − y ˉ x ˉ
∴ β 1 ^ = x y ˉ − y ˉ x ˉ x 2 ˉ − x ˉ 2 \begin{align}
\therefore \hat{\beta_1} = \frac{\bar{xy} - \bar{y}\bar{x}}{\bar{x^2} - \bar{x}^2}
\end{align} ∴ β 1 ^ = x 2 ˉ − x ˉ 2 x y ˉ − y ˉ x ˉ
と，β 1 \beta_1 β 1 β 0 \beta_0 β 0
β 1 ^ = x y ˉ − y ˉ x ˉ x 2 ˉ − x ˉ 2 = S x y S x x β 0 ^ = − β 1 ^ x ˉ + y ˉ \begin{align}
\hat{\beta_1} = \frac{\bar{xy} - \bar{y}\bar{x}}{\bar{x^2} - \bar{x}^2} = \frac{S_{xy}}{S_{xx}}\\
\hat{\beta_0} = -\hat{\beta_1}\bar{x} + \bar{y}
\end{align} β 1 ^ = x 2 ˉ − x ˉ 2 x y ˉ − y ˉ x ˉ = S xx S x y β 0 ^ = − β 1 ^ x ˉ + y ˉ
と，β 0 , β 1 \beta_0, \beta_1 β 0 , β 1
最尤法
次は最尤法での求め方です．最尤法では，誤差 ϵ i \epsilon_i ϵ i N ( 0 , σ 2 ) N(0,\sigma^2) N ( 0 , σ 2 )
すると，y i y_i y i N ( β 0 + β 1 x i , σ 2 ) N(\beta_0 + \beta_1 x_i, \sigma ^2) N ( β 0 + β 1 x i , σ 2 ) y i y_i y i
f ( y i ; β 0 , β 1 , σ 2 ) = 1 2 π σ 2 exp [ − { y i − ( β 0 + β 1 x i ) } 2 2 σ 2 ] \begin{align}
f(y_i; \beta_0, \beta_1, \sigma^2) = \frac{1}{\sqrt{2\pi\sigma^2}}\exp{[- \frac{\{y_i - (\beta_0 + \beta_1 x_i)\}^2}{2\sigma^2}]}
\end{align} f ( y i ; β 0 , β 1 , σ 2 ) = 2 π σ 2 1 exp [ − 2 σ 2 { y i − ( β 0 + β 1 x i ) } 2 ]
になります．以降，面倒なので f ( y i ; β 0 , β 1 , σ 2 ) f(y_i; \beta_0, \beta_1, \sigma^2) f ( y i ; β 0 , β 1 , σ 2 ) f ( y i ) f(y_i) f ( y i ) 尤度 ) の対数 (対数尤度関数 ) が最大になる変数 (この場合はβ 0 , β 1 \beta_0, \beta_1 β 0 , β 1
今，y i y_i y i
L ( β 0 , β 1 , σ 2 ) = ∏ i = 1 n f ( y i ) \begin{align}
L(\beta_0, \beta_1, \sigma^2) = \prod_{i=1}^n f(y_i)
\end{align} L ( β 0 , β 1 , σ 2 ) = i = 1 ∏ n f ( y i )
なので，対数尤度関数は
log ( L ( β 0 , β 1 , σ 2 ) ) = ∑ i = 1 n log f ( y i ) = − 1 2 σ 2 ∑ i = 1 n { y i − ( β 0 + β 1 x i ) } 2 − n 2 log ( 2 π σ 2 ) \begin{align}
\log(L(\beta_0, \beta_1, \sigma^2)) = \sum_{i=1}^n \log f(y_i)\\
= -\frac{1}{2\sigma^2}\sum_{i=1}^n \{y_i-(\beta_0 + \beta_1x_i) \}^2-\frac{n}{2}\log(2\pi\sigma^2)
\end{align} log ( L ( β 0 , β 1 , σ 2 )) = i = 1 ∑ n log f ( y i ) = − 2 σ 2 1 i = 1 ∑ n { y i − ( β 0 + β 1 x i ) } 2 − 2 n log ( 2 π σ 2 )
となります．あとは最小二乗法の時同様に偏微分をしていきます．すると
解法
まず対数尤度関数のうち β 0 , β 1 \beta_0, \beta_1 β 0 , β 1 β 0 , β 1 \beta_0, \beta_1 β 0 , β 1
∂ log L ∂ β 0 = 1 σ 2 ∑ i = 1 n ( y i − β 0 − β 1 x i ) = 0 \begin{align}
\frac{\partial \text{log}L}{\partial \beta_0} =\frac{1}{\sigma^2}\sum_{i=1}^n (y_i- \beta_0 - \beta_1 x_i)=0
\end{align} ∂ β 0 ∂ log L = σ 2 1 i = 1 ∑ n ( y i − β 0 − β 1 x i ) = 0
∂ log L ∂ β 1 = 1 σ 2 ∑ i = 1 n x i ( y i − β 0 − β 1 x i ) = 0 \begin{align}
\frac{\partial \text{log}L}{\partial \beta_1} =\frac{1}{\sigma^2}\sum_{i=1}^n x_i (y_i - \beta_0 - \beta_1 x_i)=0
\end{align} ∂ β 1 ∂ log L = σ 2 1 i = 1 ∑ n x i ( y i − β 0 − β 1 x i ) = 0
となり，最小二乗法の場合と同じ式になることが分かります．
次に σ 2 \sigma^2 σ 2
分かりやすくするため，一度対数尤度関数の第二項ばらします．すると
− 1 2 σ 2 ∑ i = 1 n { y i − ( β 0 + β 1 x i ) } 2 − n 2 log ( 2 π σ 2 ) = − 1 2 σ 2 ∑ i = 1 n { y i − ( β 0 + β 1 x i ) } 2 − n 2 log ( 2 π ) − n 2 log ( σ 2 ) \begin{align}
-\frac{1}{2\sigma^2}\sum_{i=1}^n \{y_i-(\beta_0 + \beta_1x_i) \}^2-\frac{n}{2}\log(2\pi\sigma^2) \\
= -\frac{1}{2\sigma^2}\sum_{i=1}^n \{y_i-(\beta_0 + \beta_1x_i) \}^2-\frac{n}{2}\log(2\pi) - \frac{n}{2}\log(\sigma^2)
\end{align} − 2 σ 2 1 i = 1 ∑ n { y i − ( β 0 + β 1 x i ) } 2 − 2 n log ( 2 π σ 2 ) = − 2 σ 2 1 i = 1 ∑ n { y i − ( β 0 + β 1 x i ) } 2 − 2 n log ( 2 π ) − 2 n log ( σ 2 )
となります．第一項は分母に σ 2 \sigma^2 σ 2 1 2 σ 4 . . . \frac{1}{2\sigma^4}... 2 σ 4 1 ... σ 2 \sigma^2 σ 2 log \log log σ 2 \sigma^2 σ 2 − n 2 σ 2 -\frac{n}{2\sigma^2} − 2 σ 2 n
よって整理すると，
∂ log L ∂ σ 2 = 1 2 σ 4 ∑ i = 1 n { y i − β 0 − x i β 1 } 2 − n 2 σ 2 = 0 \begin{align}
\frac{\partial \text{log}L}{\partial \sigma^2}= \frac{1}{2\sigma^4}\sum_{i=1}^n \{y_i - \beta_0 - x_i\beta_1\}^2 -\frac{n}{2\sigma^2}=0
\end{align} ∂ σ 2 ∂ log L = 2 σ 4 1 i = 1 ∑ n { y i − β 0 − x i β 1 } 2 − 2 σ 2 n = 0
がでてきます．最後にこれを解くと...
β 1 ^ = x y ˉ − y ˉ x ˉ x 2 ˉ − x ˉ 2 = S x y S x x , β 0 ^ = − β 1 x ˉ + y ˉ σ 2 ^ = 1 n ∑ i = 1 n ( y i − β 0 ^ − β 1 ^ x i ) 2 \begin{align}
\hat{\beta_1} = \frac{\bar{xy} - \bar{y}\bar{x}}{\bar{x^2} - \bar{x}^2} = \frac{S_{xy}}{S_{xx}}, \quad \hat{\beta_0} = -\beta_1\bar{x} + \bar{y}\\
\hat{\sigma^2} = \frac{1}{n}\sum_{i=1}^n (y_i - \hat{\beta_0}-\hat{\beta_1}x_i)^2
\end{align} β 1 ^ = x 2 ˉ − x ˉ 2 x y ˉ − y ˉ x ˉ = S xx S x y , β 0 ^ = − β 1 x ˉ + y ˉ σ 2 ^ = n 1 i = 1 ∑ n ( y i − β 0 ^ − β 1 ^ x i ) 2
が得られます．結局，回帰係数の推定値は最小二乗法の場合と一致します．
線形重回帰モデル
回帰の基本が分かったところで，モデルを多変数に拡張します．目的変数 y i , . . . , y n y_i,...,y_n y i , ... , y n x i 1 , . . . , x n p x_{i1},...,x_{np} x i 1 , ... , x n p p p p
y i = β 0 + β 1 x i 1 + . . . + β p x i p + ϵ i = β ⊤ x i + ϵ i \begin{align}
y_i &= \beta_0 + \beta_1x_{i1} + ... + \beta_{p} x_{ip} + \epsilon_i \nonumber \\
&= \bm{\beta^\top x_i} + \epsilon_i
\end{align} y i = β 0 + β 1 x i 1 + ... + β p x i p + ϵ i = β ⊤ x i + ϵ i
とベクトル表記で表せます．
さらに，y i , x i , ϵ i y_i, \bm{x_i}, \epsilon_i y i , x i , ϵ i
y = ( y 1 ⋮ y n ) , X = ( 1 x 11 ⋯ x 1 p ⋮ ⋮ ⋱ ⋮ 1 x n 1 ⋯ x n p ) , β = ( β 0 ⋮ β p ) , ϵ = ( ϵ 1 ⋮ ϵ n ) \begin{align}
\bm{y} =
\begin{pmatrix}
y_1 \\
\vdots \\
y_n\\
\end{pmatrix},
X=
\begin{pmatrix}
1 & x_{11} & \cdots & x_{1p}\\
\vdots & \vdots & \ddots & \vdots\\
1 & x_{n1} & \cdots & x_{np}\\
\end{pmatrix}, \nonumber \\
\bm{\beta} =
\begin{pmatrix}
\beta_0 \\
\vdots \\
\beta_p\\
\end{pmatrix},
\bm{\epsilon} =
\begin{pmatrix}
\epsilon_1 \\
\vdots \\
\epsilon_n\\
\end{pmatrix}\nonumber
\end{align} y = y 1 ⋮ y n , X = 1 ⋮ 1 x 11 ⋮ x n 1 ⋯ ⋱ ⋯ x 1 p ⋮ x n p , β = β 0 ⋮ β p , ϵ = ϵ 1 ⋮ ϵ n
すると先ほどの式は単に
y = X β + ϵ \begin{align}
\bm{y} = X\bm{\beta} + \bm{\epsilon}
\end{align} y = X β + ϵ
となります．やはり行列は偉大．これが線形重回帰モデル です．というかこれのことを線形回帰といったりします．何事も一般化できた方が嬉しいですもんね．この場合，β = β 0 , . . . β p \bm{\beta} = \beta_0,...\beta_p β = β 0 , ... β p 偏回帰係数 ， X X X 計画行列 と呼ばれます．
最小二乗法
ではまず最小二乗法で解いていきます．行列表記での重回帰モデルの誤差二乗和は
S ( β ) = ∥ ϵ ∥ 2 = ∥ y − X β ∥ 2 \begin{align}
S(\bm{\beta}) = \|\epsilon\|^2 = \|\bm{y} - X\bm{\beta}\|^2
\end{align} S ( β ) = ∥ ϵ ∥ 2 = ∥ y − X β ∥ 2
です．ここで ∥ ⋅ ∥ \|\cdot\| ∥ ⋅ ∥ L2ノルム (∥ x = x 1 , . . . , x n ∥ = ∑ i = 1 n x i \|\bm{x} = x_1,...,x_n\| = \sqrt{\sum_{i=1}^n x_i} ∥ x = x 1 , ... , x n ∥ = ∑ i = 1 n x i β \bm{\beta} β
∂ S ( β ) ∂ β = ∂ ∂ β ∥ y − X β ∥ 2 = 0 \begin{align}
\frac{\partial S(\bm{\beta})}{\partial \bm{\beta}} = \frac{\partial}{\partial \bm{\beta}}\|\bm{y} - X\bm{\beta}\|^2=0
\end{align} ∂ β ∂ S ( β ) = ∂ β ∂ ∥ y − X β ∥ 2 = 0
という方程式を解くと...
解法
まず，ノルムを展開します．
S ( β ) = ∥ ϵ ∥ 2 = ∥ y − X β ∥ 2 = ( y − X β ) ⊤ ( y − X β ) \begin{align}
S(\bm{\beta}) &= \|\epsilon\|^2 = \|\bm{y} - X\bm{\beta}\|^2 \nonumber\\
&= (\bm{y} - X\bm{\beta})^\top(\bm{y} - X\bm{\beta})
\end{align} S ( β ) = ∥ ϵ ∥ 2 = ∥ y − X β ∥ 2 = ( y − X β ) ⊤ ( y − X β )
整理すると
∂ S ( β ) ∂ β = ∂ ∂ β ( y − X β ) ⊤ ( y − X β ) = ∂ ∂ β [ y ⊤ y − 2 y ⊤ X β + β ⊤ X ⊤ X β ] \begin{align}
\frac{\partial S(\bm{\beta})}{\partial \bm{\beta}} &= \frac{\partial}{\partial \bm{\beta}}(\bm{y} - X\bm{\beta})^\top(\bm{y} - X\bm{\beta}) \nonumber \\
&= \frac{\partial}{\partial \bm{\beta}} [\bm{y^\top y} - 2\bm{y^\top}X\bm{\beta} + \bm{\beta^\top X^\top X \beta}]
\end{align} ∂ β ∂ S ( β ) = ∂ β ∂ ( y − X β ) ⊤ ( y − X β ) = ∂ β ∂ [ y ⊤ y − 2 y ⊤ X β + β ⊤ X ⊤ Xβ ]
ここで y ⊤ y \bm{y^\top y} y ⊤ y β \bm{\beta} β β \bm{\beta} β 行列の微分 のルールに従うと
∂ S ( β ) ∂ β = ∂ ∂ β ( − 2 y ⊤ X β + β ⊤ X ⊤ X β ) = − 2 X ⊤ y + 2 X ⊤ X β \begin{align}
\frac{\partial S(\bm{\beta})}{\partial \bm{\beta}} &=
\frac{\partial}{\partial \bm{\beta}} (- 2\bm{y^\top}X\bm{\beta} + \bm{\beta^\top X^\top X \beta}) \nonumber \\
&= -2X^\top\bm{y} + 2X^\top X \bm{\beta}
\end{align} ∂ β ∂ S ( β ) = ∂ β ∂ ( − 2 y ⊤ X β + β ⊤ X ⊤ Xβ ) = − 2 X ⊤ y + 2 X ⊤ X β
となります．最小化なので
− 2 X ⊤ y + 2 X ⊤ X β = 0 ∴ X ⊤ X β = X ⊤ y \begin{align}
-2X^\top\bm{y} + 2X^\top X \bm{\beta}=0\\
\therefore X^\top X \bm{\beta}=X^\top \bm{y}
\end{align} − 2 X ⊤ y + 2 X ⊤ X β = 0 ∴ X ⊤ X β = X ⊤ y
となります．あとはこの式を β \bm{\beta} β
進んだ注
ここでは X ⊤ X X^\top X X ⊤ X 多重共線性 が生じてはならない．また同様の理由で，n < p + 1 n < p+1 n < p + 1 オッカムの剃刀 の考え方にも通ずる．
β ^ = ( X ⊤ X ) − 1 X ⊤ y \begin{align}
\hat{\bm{\beta}} = (X^\top X)^{-1} X^{\top}\bm{y}
\end{align} β ^ = ( X ⊤ X ) − 1 X ⊤ y
と，最小二乗推定量が求まります．
最尤法
次は最尤法での解法です．単回帰の場合と同様に，ϵ ∼ N ( 0 , σ 2 ) \bm{\epsilon} \sim N(0, \sigma^2) ϵ ∼ N ( 0 , σ 2 ) y i ∼ N ( β ⊤ x i , σ 2 ) y_i \sim N(\bm{\beta^\top x_i}, \sigma^2) y i ∼ N ( β ⊤ x i , σ 2 ) y i y_i y i
f ( y i ) = 1 2 π σ 2 exp [ − ( y i − β ⊤ x i ) 2 2 σ 2 ] \begin{align}
f(y_i) = \frac{1}{\sqrt{2\pi\sigma^2}}\exp [-\frac{(y_i - \bm{\beta^\top x_i})^2}{2\sigma^2}]
\end{align} f ( y i ) = 2 π σ 2 1 exp [ − 2 σ 2 ( y i − β ⊤ x i ) 2 ]
です．対数尤度関数は
log L ( β , σ 2 ) = ∑ i = 1 n log f ( y i ) = − 1 2 σ 2 ∑ i = 1 n ( y i − β ⊤ x i ) 2 − n 2 log ( 2 π σ 2 ) = − 1 2 σ 2 ∑ i = 1 n ∥ y − X β ∥ 2 − n 2 log ( 2 π σ 2 ) \begin{align}
\log L(\bm{\beta}, \sigma^2) &= \sum_{i=1}^n \log f(y_i) \nonumber \\
&= -\frac{1}{2\sigma^2}\sum_{i=1}^n (y_i- \bm{\beta^\top x_i})^2-\frac{n}{2}\log(2\pi\sigma^2) \nonumber \\
&= -\frac{1}{2\sigma^2}\sum_{i=1}^n \|\bm{y}- X\bm{\beta}\|^2-\frac{n}{2}\log(2\pi\sigma^2)
\end{align} log L ( β , σ 2 ) = i = 1 ∑ n log f ( y i ) = − 2 σ 2 1 i = 1 ∑ n ( y i − β ⊤ x i ) 2 − 2 n log ( 2 π σ 2 ) = − 2 σ 2 1 i = 1 ∑ n ∥ y − X β ∥ 2 − 2 n log ( 2 π σ 2 )
となります．この対数尤度関数をβ \bm{\beta} β
解法
基本的にこれまでと同様です．まず尤度方程式の第二項は β \bm{\beta} β
∂ log L ( β , σ 2 ) ∂ β = − 1 2 σ 2 ∂ ∂ β ∥ y − X β ∥ 2 = − 1 2 σ 2 ∂ ∂ β ( y − X β ) ⊤ ( y − X β ) = − 1 2 σ 2 ∂ ∂ β ( y ⊤ y − 2 y ⊤ X β + β ⊤ X ⊤ X β ) \begin{align}
\frac{\partial \log L(\bm{\beta}, \sigma^2)}{\partial \bm{\beta}} &= -\frac{1}{2\sigma^2}\frac{\partial}{\partial \bm{\beta}} \|\bm{y} - X\bm{\beta}\|^2 \nonumber \\
&=-\frac{1}{2\sigma^2}\frac{\partial}{\partial \bm{\beta}} (\bm{y} - X\bm{\beta})^\top(\bm{y} - X\bm{\beta}) \nonumber \\
&= -\frac{1}{2\sigma^2}\frac{\partial}{\partial \bm{\beta}} (\bm{y^\top y} - 2\bm{y^\top}X\bm{\beta} + \bm{\beta^\top X^\top X \beta})
\end{align} ∂ β ∂ log L ( β , σ 2 ) = − 2 σ 2 1 ∂ β ∂ ∥ y − X β ∥ 2 = − 2 σ 2 1 ∂ β ∂ ( y − X β ) ⊤ ( y − X β ) = − 2 σ 2 1 ∂ β ∂ ( y ⊤ y − 2 y ⊤ X β + β ⊤ X ⊤ Xβ )
ここでやはり y ⊤ y \bm{y^\top y} y ⊤ y β \bm{\beta} β β \bm{\beta} β
∂ log L ( β , σ 2 ) ∂ β = − 1 2 σ 2 ∂ ∂ β ( − 2 y ⊤ X β + β ⊤ X ⊤ X β ) = − 1 2 σ 2 ( − 2 X ⊤ y + 2 X ⊤ X β ) = 1 σ 2 X ⊤ ( y − X β ) \begin{align}
\frac{\partial \log L(\bm{\beta}, \sigma^2)}{\partial \bm{\beta}} &=
-\frac{1}{2\sigma^2} \frac{\partial}{\partial \bm{\beta}} (- 2\bm{y^\top}X\bm{\beta} + \bm{\beta^\top X^\top X \beta}) \nonumber \\
&= -\frac{1}{2\sigma^2}(-2X^\top\bm{y} + 2X^\top X \bm{\beta}) \nonumber \\
&= \frac{1}{\sigma^2}X^\top (\bm{y} - X\bm{\beta})
\end{align} ∂ β ∂ log L ( β , σ 2 ) = − 2 σ 2 1 ∂ β ∂ ( − 2 y ⊤ X β + β ⊤ X ⊤ Xβ ) = − 2 σ 2 1 ( − 2 X ⊤ y + 2 X ⊤ X β ) = σ 2 1 X ⊤ ( y − X β )
となります．最小化なので
1 σ 2 X ⊤ ( y − X β ) = 0 ∴ X ⊤ X β = X ⊤ y \begin{align}
\frac{1}{\sigma^2}X^\top (\bm{y} - X\bm{\beta})=0\\
\therefore X^\top X \bm{\beta}=X^\top \bm{y}
\end{align} σ 2 1 X ⊤ ( y − X β ) = 0 ∴ X ⊤ X β = X ⊤ y
となります．あとはこの式を β \bm{\beta} β
β ^ = ( X ⊤ X ) − 1 X ⊤ y \begin{align}
\hat{\bm{\beta}} = (X^\top X)^{-1} X^{\top}\bm{y}
\end{align} β ^ = ( X ⊤ X ) − 1 X ⊤ y
が求まります．次に， σ 2 \sigma^2 σ 2
log L ( β , σ 2 ) = ∑ i = 1 n log f ( y i ) = − 1 2 σ 2 ∑ i = 1 n ∥ y − X β ∥ 2 − n 2 log ( 2 π σ 2 ) \begin{align}
\log L(\bm{\beta}, \sigma^2) &= \sum_{i=1}^n \log f(y_i) \nonumber \\
&= -\frac{1}{2\sigma^2}\sum_{i=1}^n \|\bm{y}- X\bm{\beta}\|^2-\frac{n}{2}\log(2\pi\sigma^2)
\end{align} log L ( β , σ 2 ) = i = 1 ∑ n log f ( y i ) = − 2 σ 2 1 i = 1 ∑ n ∥ y − X β ∥ 2 − 2 n log ( 2 π σ 2 )
の第一，第二どちらの項もσ 2 \sigma^2 σ 2
∂ log L ( β , σ 2 ) ∂ β = 1 2 σ 4 ∑ i = 1 n ∥ y − X β ∥ 2 − n 2 σ 2 \begin{align}
\frac{\partial \log L(\bm{\beta}, \sigma^2)}{\partial \bm{\beta}} = \frac{1}{2\sigma^4}\sum_{i=1}^n \|\bm{y}- X\bm{\beta}\|^2-\frac{n}{2\sigma^2}
\end{align} ∂ β ∂ log L ( β , σ 2 ) = 2 σ 4 1 i = 1 ∑ n ∥ y − X β ∥ 2 − 2 σ 2 n
となります．これを0として解くと
σ 2 ^ = 1 n ∥ y − X β ^ ∥ 2 \begin{align}
\hat{\sigma^2} = \frac{1}{n}\| \bm{y} - X \hat{\bm{\beta}}\|^2
\end{align} σ 2 ^ = n 1 ∥ y − X β ^ ∥ 2
となります．よってまとめると，
β ^ = ( X ⊤ X ) − 1 X ⊤ y , σ 2 ^ = 1 n ∥ y − X β ^ ∥ 2 \begin{align}
\hat{\bm{\beta}} = (X^\top X)^{-1} X^{\top}\bm{y}, \quad
\hat{\sigma^2} = \frac{1}{n}\| \bm{y} - X \hat{\bm{\beta}}\|^2
\end{align} β ^ = ( X ⊤ X ) − 1 X ⊤ y , σ 2 ^ = n 1 ∥ y − X β ^ ∥ 2
と，β \bm{\beta} β σ 2 \sigma^2 σ 2
多項式回帰
多項式回帰は，目的変数 y i y_i y i x i x_i x i n n n
これまで見てきた線形回帰は，どれもデータの分布が線形なことを前提としていました．多項式回帰を含む非線形のモデルでは，より曲線的なデータ分布に対しても上手くあてはまるモデルを考えることが出来ます．
まずは例を見てみます．
pythonコード
np. random. seed( 0 )
X = 2 - 3 * np. random. normal( 0 , 1 , 100 )
y = X - 2 * ( X ** 2 ) + 0.5 * ( X ** 3 ) + np. random. normal( - 3 , 3 , 100 )
X = X[ : , np. newaxis]
plt. scatter( X, y, color= 'blue' , label= 'Data Points' )
degrees = [ 1 , 2 , 3 ]
colors = [ 'red' , 'green' , 'orange' ]
labels = [ 'Degree 1' , 'Degree 2' , 'Degree 3' ]
X_plot = np. linspace( X. min ( ) , X. max ( ) , 100 )
X_plot = X_plot[ : , np. newaxis]
for i, degree in enumerate ( degrees) :
polynomial_features = PolynomialFeatures( degree= degree)
X_poly = polynomial_features. fit_transform( X)
model = LinearRegression( )
model. fit( X_poly, y)
y_poly_plot = model. predict( polynomial_features. fit_transform( X_plot) )
plt. plot( X_plot, y_poly_plot, color= colors[ i] , label= labels[ i] )
plt. xlabel( 'X' , fontsize= 15 )
plt. ylabel( 'y' , fontsize= 15 )
plt. legend( )
plt. grid( )
plt. title( 'Polynomial Regression with Different Degrees' , fontsize= 20 )
plt. savefig( '../figures/polynomial.png' )
この図の青点がデータの分布です．これに対し，赤緑橙の曲線が多項式回帰の結果です．橙はよく当てはまっている一方，赤や緑は微妙です．
この差は，回帰に使う項の数となります．重回帰の場合には，x i x_i x i x n x_n x n
今は 1 1 1 x x x
多項式回帰のモデルは
y = β 0 + β 1 x + β 2 x 2 + ⋯ β n x n + ϵ \begin{align}
y = \beta_0 + \beta_1 x + \beta_2 x^2 + \cdots \beta_n x^n + \epsilon
\end{align} y = β 0 + β 1 x + β 2 x 2 + ⋯ β n x n + ϵ
です．式の形は重回帰の場合と同じなので，同様にして最小二乗推定値
β ^ = ( X ⊤ X ) − 1 X ⊤ y \begin{align}
\hat{\bm{\beta}} = (X^\top X)^{-1} X^{\top}\bm{y}
\end{align} β ^ = ( X ⊤ X ) − 1 X ⊤ y
を求めることができます．また，これまた重回帰の場合同様に項の数は多ければよいというものでもありません．適切な値を選ぶ必要があります．
