線形回帰のベクトル・行列による表現

設定

以下のようなp個の変数 (x_1, x_2, ..., x_{p-1}) による線形回帰を考える：

\begin{aligned} y_i &= \beta_0 + \beta_1 x_{1i} + \beta_2 x_{2i} + ... \beta_p x_{p-1,i} + \varepsilon_i\\ &= \boldsymbol{x}_i^\top \boldsymbol{\beta} + \varepsilon_i. \end{aligned}

ここで、

\boldsymbol{x}_i = \begin{pmatrix} 1\\ x_{1i}\\ x_{2i}\\ \vdots \\ x_{p-1, i}\\ \end{pmatrix}, \quad \boldsymbol{\beta} = \begin{pmatrix} \beta_0\\ \beta_1\\ \beta_2\\ \vdots \\ \beta_{p-1}\\ \end{pmatrix}

のようにおいた。
なお、誤差項\varepsilon_i は以下のように分散が\sigma^2 の独立同時なホワイトノイズである：

E[\varepsilon_i]=0, \quad E[\varepsilon_{i}\varepsilon_{j}] = \left\{ \begin{array}{ll} \sigma^2 & i=j \\ 0 & i \ne j \end{array} \right. .

観測値の行列を用いた表現

N個の観測値\{ (\boldsymbol{x}_1, y_1), (\boldsymbol{x}_2, y_2), ..., (\boldsymbol{x}_N, y_N) \} が与えられているとき、

\boldsymbol{Y} = \begin{pmatrix} y_1\\ y_2\\ \vdots \\ y_N \end{pmatrix}, \quad \boldsymbol{\varepsilon} = \begin{pmatrix} \varepsilon_1\\ \varepsilon_2\\ \vdots \\ \varepsilon_N \end{pmatrix},

\boldsymbol{X} = \begin{pmatrix} \boldsymbol{x}^\top_1\\ \boldsymbol{x}^\top_2\\ \vdots \\ \boldsymbol{x}^\top_N \end{pmatrix} = \begin{pmatrix} 1 & x_{11} & x_{21} & ... & x_{p-1,1}\\ 1 & x_{12} & x_{22} & ... & x_{p-1,2}\\ \vdots & & & \vdots & \vdots\\ 1 & x_{1N} & x_{2N} & ... & x_{p-1,N} \end{pmatrix}

とすると、回帰式は以下のように表せる：

\boldsymbol{Y} = \boldsymbol{X} \boldsymbol{\beta} + \boldsymbol{\varepsilon}

例： p=2（一変数、単回帰）の場合

y_i = \beta_0 + \beta_1 x_i + \varepsilon_i

の場合を考える。このとき、

\boldsymbol{X} = \begin{pmatrix} 1 & x_1 \\ 1 & x_2\\ \vdots \\ 1 & x_N \end{pmatrix}, \quad \boldsymbol{\beta} = \begin{pmatrix} \beta_0\\ \beta_1 \end{pmatrix}

であるので、\boldsymbol{Y} = \boldsymbol{X} \boldsymbol{\beta} + \boldsymbol{\varepsilon} を成分ごとに表すと、

\begin{aligned} \begin{pmatrix} y_1\\ y_2\\ \vdots \\ y_N \end{pmatrix} &= \begin{pmatrix} 1 & x_1 \\ 1 & x_2\\ \vdots \\ 1 & x_N \end{pmatrix} \begin{pmatrix} \beta_0\\ \beta_1 \end{pmatrix} + \begin{pmatrix} \varepsilon_1\\ \varepsilon_2\\ \vdots \\ \varepsilon_N \end{pmatrix}\\ &= \begin{pmatrix} \beta_0 + \beta_1 x_1 + \varepsilon_1\\ \beta_0 + \beta_1 x_2 + \varepsilon_2\\ \vdots \\ \beta_0 + \beta_1 x_N + \varepsilon_N\\ \end{pmatrix} \end{aligned}

線形回帰における係数推定量のベクトル・行列による表現

回帰式

\boldsymbol{Y} = \boldsymbol{X} \boldsymbol{\beta} + \boldsymbol{\varepsilon}

が与えられた際の係数 \boldsymbol{\beta} の最小二乗推定量は、以下のように表される：

\begin{aligned} \hat{\boldsymbol{\beta}} &= (\boldsymbol{X}^\top \boldsymbol{X})^{-1} \boldsymbol{X}^\top \boldsymbol{Y} \\ &= \left(\sum_{i=1}^N \boldsymbol{x}_i \boldsymbol{x}_i^\top\right)^{-1} \left( \sum_{i=1}^N \boldsymbol{x}_i y_i \right) \tag{1} \end{aligned}

この時、説明変数\boldsymbol{X}と誤差項\boldsymbol{\varepsilon}が無相関であれば、この推定量は不偏推定量になる：

E[\hat{\boldsymbol{\beta}}] = \boldsymbol{\beta}. \tag{2}

また、このとき推定量\hat{\boldsymbol{\beta}}の分散共分散行列は、以下のように表される：

\begin{aligned} V[\hat{\boldsymbol{\beta}}] &= \sigma^2 (\boldsymbol{X}^\top \boldsymbol{X})^{-1}\\ &= \sigma^2 \left(\sum_{i=1}^N \boldsymbol{x}_i \boldsymbol{x}_i^\top\right)^{-1}. \tag{3} \end{aligned}

特に誤差項 \varepsilon_t が正規分布 \varepsilon \sim N(0, \sigma^2) に従う場合、
\hat{\boldsymbol{\beta}} の分布も正規分布になる：

\hat{\boldsymbol{\beta}} \sim N(\boldsymbol{\beta}, \sigma^2 (\boldsymbol{X}^\top \boldsymbol{X})^{-1})

例：p=2（一変数、単回帰）の場合

y_i = \beta_0 + \beta_1 x_i + \varepsilon_i

の場合を考える。

係数推定量

このとき、

\boldsymbol{X} = \begin{pmatrix} 1 & x_1 \\ 1 & x_2\\ \vdots \\ 1 & x_N \end{pmatrix}, \quad \boldsymbol{\beta} = \begin{pmatrix} \beta_0\\ \beta_1 \end{pmatrix}

であるので、

\begin{aligned} (\boldsymbol{X}^\top \boldsymbol{X})^{-1} &= \begin{pmatrix} \sum_i 1 & \sum_i x_i \\ \sum_i x_i & \sum_i x_i^2 \end{pmatrix}^{-1}\\ &=\frac{1}{N \sum_i x_i^2 - \left( \sum_i x_i \right)^2} \begin{pmatrix} \sum_i x_i^2 & -\sum_i x_i \\ -\sum_i x_i & N \end{pmatrix}\\ &= \frac{1}{\sum_i (x_i - \bar{x})^2} \begin{pmatrix} \frac{1}{N}\sum_i x_i^2 & - \bar{x} \\ - \bar{x} & 1 \end{pmatrix} \end{aligned}

これを\hat{\boldsymbol{\beta}} = (\boldsymbol{X}^\top \boldsymbol{X})^{-1} \boldsymbol{X}^\top \boldsymbol{Y}に代入することで、

\begin{aligned} \hat{\boldsymbol{\beta}} &= \frac{1}{\sum_i (x_i - \bar{x})^2} \begin{pmatrix} \frac{1}{N}\sum_i x_i^2 & - \bar{x} \\ - \bar{x} & 1 \end{pmatrix} \begin{pmatrix} \sum_i y_i \\ \sum_i x_i y_i \end{pmatrix}\\ \end{aligned}

が得られる（\sum_iは\sum_{i=1}^Nの略記、\bar{x} = \frac{1}{N} \sum_{i=1}^N x_i）。

これを整理すると、

\begin{aligned} \hat{\beta}_0 &= \frac{\bar{y} \sum_{i=1}^N x_i^2 - \bar{x} \sum_{i=1}^N x_i y_i }{\sum_{i=1}^N (x_i - \bar{x})^2}\\ &= \bar{y} - \hat{\beta}_1 \bar{x}\\ \hat{\beta}_1 &= \frac{\sum_{i=1}^N (y_i - \bar{y})(x_i - \bar{x})}{\sum_{i=1}^N (x_i - \bar{x})^2} \end{aligned}

ただし、\bar{y}=\frac{1}{N} \sum_{i=1}^N y_iとした。

係数推定量の分散・共分散

前述のように

(\boldsymbol{X}^\top \boldsymbol{X})^{-1} = \frac{1}{\sum_i (x_i - \bar{x})^2} \begin{pmatrix} \frac{1}{N}\sum_i x_i^2 & - \bar{x} \\ - \bar{x} & 1 \end{pmatrix}

であったので、これを V[\hat{\boldsymbol{\beta}}] = \sigma^2 (\boldsymbol{X}^\top \boldsymbol{X})^{-1} に代入すると、

\begin{aligned} V[\hat{\beta}_0] &= \frac{\sigma^2 \frac{1}{N} \sum_{i=1}^N x_i^2}{\sum_{i=1}^N (x_i - \bar{x})^2}\\ V[\hat{\beta}_1] &= \frac{\sigma^2 }{\sum_{i=1}^N (x_i - \bar{x})^2}\\ Cov[\hat{\beta}_0, \hat{\beta}_1] &= -\frac{\sigma^2 \bar{x}}{\sum_{i=1}^N (x_i - \bar{x})^2} \end{aligned}

係数推定量の分布

以上を\hat{\boldsymbol{\beta}} \sim N(\boldsymbol{\beta}, \sigma^2 (\boldsymbol{X}^\top \boldsymbol{X})^{-1})に代入して、

\begin{pmatrix} \hat{c}\\ \hat{\phi} \end{pmatrix} \sim N \left( \begin{pmatrix} c\\ \phi \end{pmatrix}, \frac{\sigma^2}{\sum_i (x_i - \bar{x})^2} \begin{pmatrix} \frac{1}{N}\sum_i x_i^2 & - \bar{x} \\ - \bar{x} & 1 \end{pmatrix} \right)

導出・証明

(1) の導出

\hat{\boldsymbol{\beta}} = (\boldsymbol{X}^\top \boldsymbol{X})^{-1} \boldsymbol{X}^\top \boldsymbol{Y} を示す。

残差二乗和は以下のように表される：

RSS(\boldsymbol{\beta}) = (\boldsymbol{Y}-\boldsymbol{X}\boldsymbol{\beta})^\top (\boldsymbol{Y}-\boldsymbol{X}\boldsymbol{\beta})

これを \boldsymbol{\beta} で微分すると、

\begin{aligned} \frac{\partial RSS(\boldsymbol{\beta})}{\partial \boldsymbol{\beta}} &= \frac{\partial (\boldsymbol{Y}^\top \boldsymbol{Y} - \boldsymbol{Y}^\top \boldsymbol{X}\boldsymbol{\beta} - \boldsymbol{\beta}^\top \boldsymbol{X}^\top \boldsymbol{Y} + \boldsymbol{\beta}^\top \boldsymbol{X}^\top \boldsymbol{X} \boldsymbol{\beta})}{\partial \boldsymbol{\beta}} \\ &= \frac{\partial (\boldsymbol{Y}^\top \boldsymbol{Y})}{\partial \boldsymbol{\beta}} - 2 \frac{\partial (\boldsymbol{Y}^\top \boldsymbol{X}\boldsymbol{\beta})}{\partial \boldsymbol{\beta}} + \frac{\partial (\boldsymbol{\beta}^\top \boldsymbol{X}^\top \boldsymbol{X} \boldsymbol{\beta})}{\partial \boldsymbol{\beta}}\\ &= \boldsymbol{0} - 2 (\boldsymbol{Y}^\top \boldsymbol{X})^\top + \{ \boldsymbol{X}^\top \boldsymbol{X} + (\boldsymbol{X}^\top \boldsymbol{X})^\top \} \boldsymbol{\beta}\\ &= -2 \boldsymbol{X}^\top \boldsymbol{Y} + 2 \boldsymbol{X}^\top \boldsymbol{X} \boldsymbol{\beta} \end{aligned}

のようになる。
1行目から2行目では、\boldsymbol{Y}^\top \boldsymbol{X}\boldsymbol{\beta} がスカラーであることから\boldsymbol{\beta}^\top \boldsymbol{X}^\top \boldsymbol{Y} = (\boldsymbol{Y}^\top \boldsymbol{X}\boldsymbol{\beta})^\top = \boldsymbol{Y}^\top \boldsymbol{X}\boldsymbol{\beta}を用いた。
2行目から3行目では、以下のベクトルによる微分の関係式を適用した：

\frac{\partial (\boldsymbol{A}\boldsymbol{x})}{\partial \boldsymbol{x}} = \boldsymbol{A}^\top, \quad \frac{\partial (\boldsymbol{x}^\top \boldsymbol{A}\boldsymbol{x})}{\partial \boldsymbol{x}} = (\boldsymbol{A} + \boldsymbol{A}^\top) \boldsymbol{x}.

RSS(\boldsymbol{\beta}) が最小になるとき上記の微分の値は\boldsymbol{0} になるので、

\frac{\partial RSS(\boldsymbol{\beta})}{\partial \boldsymbol{\beta}} \Bigr|_{\boldsymbol{\beta}=\hat{\boldsymbol{\beta}}} = \boldsymbol{0}

より、

\boldsymbol{X}^\top \boldsymbol{X} \hat{\boldsymbol{\beta}} = \boldsymbol{X}^\top \boldsymbol{Y}

を解いて、

\hat{\boldsymbol{\beta}} = (\boldsymbol{X}^\top \boldsymbol{X})^{-1} \boldsymbol{X}^\top \boldsymbol{Y}

(2) の証明

E[\hat{\boldsymbol{\beta}}] = \boldsymbol{\beta} を示す。

\begin{aligned} \hat{\boldsymbol{\beta}} &= (\boldsymbol{X}^\top \boldsymbol{X})^{-1} \boldsymbol{X}^\top \boldsymbol{Y}\\ &= (\boldsymbol{X}^\top \boldsymbol{X})^{-1} \boldsymbol{X}^\top (\boldsymbol{X} \boldsymbol{\beta} + \boldsymbol{\varepsilon})\\ &= (\boldsymbol{X}^\top \boldsymbol{X})^{-1} \boldsymbol{X}^\top \boldsymbol{X} \boldsymbol{\beta} + (\boldsymbol{X}^\top \boldsymbol{X})^{-1} \boldsymbol{X}^\top \boldsymbol{\varepsilon}\\ &= \boldsymbol{\beta} + (\boldsymbol{X}^\top \boldsymbol{X})^{-1} \boldsymbol{X}^\top \boldsymbol{\varepsilon} \\ &(\because (\boldsymbol{X}^\top \boldsymbol{X})^{-1} (\boldsymbol{X}^\top \boldsymbol{X} ) = \boldsymbol{1}) \end{aligned}

ここで、\boldsymbol{X}と\boldsymbol{\varepsilon}が無相関より、

\begin{aligned} E[\hat{\boldsymbol{\beta}}] &= \boldsymbol{\beta} + (\boldsymbol{X}^\top \boldsymbol{X})^{-1} \boldsymbol{X}^\top E[ \boldsymbol{\varepsilon}]\\ &= \boldsymbol{\beta} + \boldsymbol{0}\\ &= \boldsymbol{\beta} \end{aligned}

(3)の証明

V[\hat{\boldsymbol{\beta}}] = \sigma^2 (\boldsymbol{X}^\top \boldsymbol{X})^{-1} を示す。

\begin{aligned} V[\hat{\boldsymbol{\beta}}] &= E[( \hat{\boldsymbol{\beta}} - \boldsymbol{\beta})( \hat{\boldsymbol{\beta}} - \boldsymbol{\beta})^\top]\\ &= E[\{(\boldsymbol{X}^\top \boldsymbol{X})^{-1} \boldsymbol{X}^\top \boldsymbol{\varepsilon} \}\{(\boldsymbol{X}^\top \boldsymbol{X})^{-1} \boldsymbol{X}^\top \boldsymbol{\varepsilon} \}^\top]\\ &= (\boldsymbol{X}^\top \boldsymbol{X})^{-1} \boldsymbol{X}^\top E[\boldsymbol{\varepsilon}\boldsymbol{\varepsilon}^\top] (\boldsymbol{X}^\top)^\top \{(\boldsymbol{X}^\top \boldsymbol{X})^{-1}\}^\top\\ &= (\boldsymbol{X}^\top \boldsymbol{X})^{-1} \boldsymbol{X}^\top \sigma^2 \boldsymbol{1} \boldsymbol{X} (\boldsymbol{X}^\top \boldsymbol{X})^{-1}\\ &= \sigma^2 (\boldsymbol{X}^\top \boldsymbol{X})^{-1} (\boldsymbol{X}^\top \boldsymbol{X} )(\boldsymbol{X}^\top \boldsymbol{X})^{-1}\\ &= \sigma^2 (\boldsymbol{X}^\top \boldsymbol{X})^{-1} \end{aligned}

補足

一般のpで\boldsymbol{X}^\top \boldsymbol{X}を成分表示する。

\boldsymbol{x}_i = \begin{pmatrix} 1\\ x_{1i}\\ x_{2i}\\ \vdots \\ x_{pi}\\ \end{pmatrix} , \quad \boldsymbol{X} = \begin{pmatrix} \boldsymbol{x}^\top_1\\ \boldsymbol{x}^\top_2\\ \vdots \\ \boldsymbol{x}^\top_N \end{pmatrix} = \begin{pmatrix} 1 & x_{11} & x_{21} & ... & x_{p1}\\ 1 & x_{12} & x_{22} & ... & x_{p2}\\ \vdots & & & \vdots & \vdots\\ 1 & x_{1N} & x_{2N} & ... & x_{pN} \end{pmatrix}

であるので、

\begin{aligned} \boldsymbol{X}^\top \boldsymbol{X} &= \sum_{i=1}^N \boldsymbol{x}_i \boldsymbol{x}_i^\top\\ &= \begin{pmatrix} \sum_i 1 & \sum_i x_{1i} & \sum_i x_{2i} & ... & \sum_i x_{pi}\\ \sum_i x_{1i} & \sum_i x_{1i}^2 & \sum_i x_{1i} x_{2i} & ... & \sum_i x_{1i} x_{pi}\\ \sum_i x_{2i} & \sum_i x_{2i}x_{1i} & \sum_i x_{2i}^2 & ... & \sum_i x_{2i} x_{pi}\\ \vdots & & & \vdots & \vdots\\ \sum_i x_{pi} & \sum_i x_{pi}x_{1i} & \sum_i x_{pi}x_{2i} & ... & \sum_i x_{pi}^2 & \end{pmatrix} \end{aligned}

のように表される。

Ridge 回帰とその推定量の性質

以下の投稿で書いたように、線形回帰 \boldsymbol{Y} = \boldsymbol{X}\boldsymbol{\beta} + \boldsymbol{\varepsilon} の回帰係数 \boldsymbol{\beta} の推定量 \hat{\boldsymbol{\beta}} は、

\hat{\boldsymbol{\beta}} = (\boldsymbol{X}^\top \boldsymbol{X})^{-1} \boldsymbol{X}^\top \boldsymbol{Y} \tag{1}

のように表すことができた。

https://zenn.dev/link/comments/db9e98ef7a0301

しかし、上記の式が成立するためには、 \boldsymbol{X}^\top \boldsymbol{X} が逆行列を持たなくてはいけない。

一方、Ridge 回帰では、\boldsymbol{X}^\top \boldsymbol{X} が逆行列を持たないような場合でも常に係数 \boldsymbol{\beta} の推定を行うことができる。

本投稿では Ridge 回帰とその推定量の性質について解説する。

結果

サンプルサイズを n、説明変数の数を p とする。
このとき、\boldsymbol{X} \in \mathbb{R}^{n \times p}, \, \boldsymbol{Y} \in \mathbb{R}^{n \times 1}, \, \boldsymbol{\beta} \in \mathbb{R}^{p \times 1} である。

Ridge 回帰による 係数 \boldsymbol{\beta} の推定量は、以下である：

\hat{\boldsymbol{\beta}}_\alpha = (\boldsymbol{X}^\top \boldsymbol{X} + \alpha I_p)^{-1} \boldsymbol{X}^\top \boldsymbol{Y} \tag{2}

I_p は p \times p の単位行列である。

この推定量は、 \alpha > 0 である限り常に存在する。

なお、\hat{\boldsymbol{\beta}}_\alpha は \boldsymbol{\beta} の不偏推定量にはならない（E[\hat{\boldsymbol{\beta}}_\alpha] \ne \boldsymbol{\beta}）。

Ridge 回帰の概要

Ridge 回帰では、推定量を求める際に残差二乗和 ||\boldsymbol{Y} - \boldsymbol{X} \boldsymbol{\beta}||^2 そのものではなく係数 \boldsymbol{\beta} の二乗ノルムを足したもの

L(\boldsymbol{\beta}) = ||\boldsymbol{Y} - \boldsymbol{X} \boldsymbol{\beta}||^2 + \alpha ||\boldsymbol{\beta}||^2 \quad (\alpha>0: {\rm const.})

を最小化する。
こうすることで、\boldsymbol{X}^\top \boldsymbol{X} が逆行列を持たない場合でも推定を行えるほか、係数 \boldsymbol{\beta} の推定量が不安定さに起因して爆発的に大きくなることも抑えられる。

Ridge 回帰による推定量の導出

上記の正則化項付きの損失関数L(\boldsymbol{\beta})を最小化することで、冒頭の式(2)で示した推定量 \hat{\boldsymbol{\beta}}_\alpha が求まる。

これを示すために、L(\boldsymbol{\beta}) を \boldsymbol{\beta} で微分する：

\begin{aligned} \frac{\partial L(\boldsymbol{\beta})}{\partial \boldsymbol{\beta}} &= - 2 \boldsymbol{X}^\top (\boldsymbol{Y} - \boldsymbol{X} \boldsymbol{\beta} ) + 2 \alpha \boldsymbol{\beta}\\ &= -2 \boldsymbol{X}^\top \boldsymbol{Y} + 2 (\boldsymbol{X}^\top \boldsymbol{X} + \alpha I_p)\boldsymbol{\beta}. \end{aligned}

L(\boldsymbol{\beta}) は最小になる\boldsymbol{\beta} = \hat{\boldsymbol{\beta}}_\alpha では、上記の微分の値が \boldsymbol{0} になる。

従ってこのとき、

(\boldsymbol{X}^\top \boldsymbol{X} + \alpha I_p) \hat{\boldsymbol{\beta}}_\alpha = \boldsymbol{X}^\top \boldsymbol{Y}

が成り立つ。

ここで、行列 (\boldsymbol{X}^\top \boldsymbol{X} + \alpha I_p) は \alpha > 0 で常に非特異で逆行列を持つ^[1]ので、先ほどの式の両辺に左から (\boldsymbol{X}^\top \boldsymbol{X} + \alpha I_p)^{-1} をかけることで式(2)が求まる。

Ridge 回帰推定量の期待値

通常の最小二乗推定では、推定量 \hat{\boldsymbol{\beta}} の期待値は真の値 \boldsymbol{\beta} と一致することから、\hat{\boldsymbol{\beta}} は不偏推定量であった。

しかし、Ridge 回帰による推定量 \hat{\boldsymbol{\beta}}_\alpha は不偏推定量にならない（E[\hat{\boldsymbol{\beta}}_\alpha] \ne \boldsymbol{\beta}）ことを示す。

まず、式(2) の右辺に \boldsymbol{Y} = \boldsymbol{X}\boldsymbol{\beta} + \boldsymbol{\varepsilon} を代入する：

\begin{aligned} \hat{\boldsymbol{\beta}}_\alpha &= (\boldsymbol{X}^\top \boldsymbol{X} + \alpha I_p)^{-1} \boldsymbol{X}^\top (\boldsymbol{X}\boldsymbol{\beta} + \boldsymbol{\varepsilon})\\ &= (\boldsymbol{X}^\top \boldsymbol{X} + \alpha I_p)^{-1} \boldsymbol{X}^\top \boldsymbol{X} \boldsymbol{\beta} + (\boldsymbol{X}^\top \boldsymbol{X} + \alpha I_p)^{-1} \boldsymbol{X}^\top \boldsymbol{\varepsilon} \tag{3} \end{aligned}

この両辺について期待値をとると、

\begin{aligned} E[\hat{\boldsymbol{\beta}}_\alpha] &= (\boldsymbol{X}^\top \boldsymbol{X} + \alpha I_p)^{-1} \boldsymbol{X}^\top \boldsymbol{X} \boldsymbol{\beta}\\ &= (\boldsymbol{X}^\top \boldsymbol{X} + \alpha I_p)^{-1} (\boldsymbol{X}^\top \boldsymbol{X} + \alpha I_p) \boldsymbol{\beta} - (\boldsymbol{X}^\top \boldsymbol{X} + \alpha I_p)^{-1} \alpha I_p \boldsymbol{\beta}\\ &= \boldsymbol{\beta} - \alpha (\boldsymbol{X}^\top \boldsymbol{X} + \alpha I_p)^{-1} \boldsymbol{\beta}. \end{aligned}

以上から、\hat{\boldsymbol{\beta}}_\alpha の期待値は \boldsymbol{\beta} と一致せず、- \alpha (\boldsymbol{X}^\top \boldsymbol{X} + \alpha I_p)^{-1} \boldsymbol{\beta} だけバイアスが入る。

参考： Ridge 回帰推定量の分散

参考までに、Ridge 回帰推定量 \hat{\boldsymbol{\beta}}_\alpha の分散についても見てみる。

式(3) より、

\begin{aligned} V[\hat{\boldsymbol{\beta}}_\alpha] &= E[\{ (\boldsymbol{X}^\top \boldsymbol{X} + \alpha I_p)^{-1} \boldsymbol{X}^\top \boldsymbol{\varepsilon} \} \{ (\boldsymbol{X}^\top \boldsymbol{X} + \alpha I_p)^{-1} \boldsymbol{X}^\top \boldsymbol{\varepsilon}\}^\top]\\ &= (\boldsymbol{X}^\top \boldsymbol{X} + \alpha I_p)^{-1} \boldsymbol{X}^\top E[\boldsymbol{\varepsilon} \boldsymbol{\varepsilon}^\top] \boldsymbol{X} \{ (\boldsymbol{X}^\top \boldsymbol{X} + \alpha I_p)^{-1} \}^\top\\ &= \sigma^2 (\boldsymbol{X}^\top \boldsymbol{X} + \alpha I_p)^{-1} \boldsymbol{X}^\top \boldsymbol{X} (\boldsymbol{X}^\top \boldsymbol{X} + \alpha I_p)^{-1} \end{aligned}

となる。

2行目から3行目での変形では、

E[\boldsymbol{\varepsilon} \boldsymbol{\varepsilon}^\top] = \sigma^2 I_n

\{ (\boldsymbol{X}^\top \boldsymbol{X} + \alpha I_p)^{-1} \}^\top = (\boldsymbol{X}^\top \boldsymbol{X} + \alpha I_p)^{-1}

であることを用いた。
後者については、(\boldsymbol{X}^\top \boldsymbol{X} + \alpha I_p) が対称行列であることと、対称行列の逆行列もまた対称行列であることに着目した。

脚注

1. https://zenn.dev/link/comments/e9090f5ef1defc の補足参照 ↩︎

残差二乗和の統計的性質と誤差項分散

以下のような線形回帰モデルを考える：

\boldsymbol{Y} = \boldsymbol{X}\boldsymbol{\beta} + \boldsymbol{\varepsilon}.

ここで、\boldsymbol{Y} \in \mathbb{R}^n は目的変数ベクトル、\boldsymbol{X} \in \mathbb{R}^{n\times p} は p 個の説明変数から成る計画行列、\boldsymbol{\beta} \in \mathbb{R}^p は回帰係数、\boldsymbol{\varepsilon} \in \mathbb{R}^n は誤差項である。
誤差項 \boldsymbol{\varepsilon} = (\varepsilon_1, \varepsilon_2, ... \varepsilon_n)^\top は独立同一な分布に従いその期待値は0、分散は \sigma^2 とする。

計画行列 \boldsymbol{X} は full rank であるものとして、回帰係数の最小二乗法による推定量は以下のように表すことができる：

\hat{\boldsymbol{\beta}} = (\boldsymbol{X}^\top \boldsymbol{X})^{-1} \boldsymbol{X}^\top \boldsymbol{Y}.

これを用いて、残差 \hat{\boldsymbol{\varepsilon}} を以下のように定義する：

\begin{aligned} \hat{\boldsymbol{\varepsilon}} &= \boldsymbol{Y} - \boldsymbol{X} \hat{\boldsymbol{\beta}}\\ &= \left\{\boldsymbol{I}_n - \boldsymbol{X} (\boldsymbol{X}^\top \boldsymbol{X})^{-1} \boldsymbol{X}^\top \right\}\boldsymbol{Y}. \end{aligned}

以降、 \boldsymbol{M}_X = \boldsymbol{I}_n - \boldsymbol{X} (\boldsymbol{X}^\top \boldsymbol{X})^{-1} \boldsymbol{X}^\top と置く。

残差二乗和 \hat{\boldsymbol{\varepsilon}}^\top\hat{\boldsymbol{\varepsilon}} の統計的性質

残差二乗和 \hat{\boldsymbol{\varepsilon}}^\top\hat{\boldsymbol{\varepsilon}} について、

s^2 = \frac{\hat{\boldsymbol{\varepsilon}}^\top\hat{\boldsymbol{\varepsilon}}}{n-p}

とおくと、s^2 は誤差項分散 \sigma^2 の不偏推定量になる：

E[s^2] = \sigma^2

特に誤差項 \{ \varepsilon_i \}_{i=1}^n が正規分布 \varepsilon_i \sim N(0, \sigma^2) に従うとき、

\frac{\hat{\boldsymbol{\varepsilon}}^\top\hat{\boldsymbol{\varepsilon}}}{\sigma^2} \sim \chi^2(n-p)

が成り立つ。
カイ二乗分布 \chi^2(k) の平均値 k であるので E\left[\hat{\boldsymbol{\varepsilon}}^\top\hat{\boldsymbol{\varepsilon}} / \sigma^2 \right] = n-p であることがわかり、このことからも E[s^2] = \sigma^2 であることがわかる。

なお、\chi^2(k) の分散は 2k であるため、推定量 s^2 の分散は以下のように表すことができる：

\begin{aligned} V[s^2] &= \frac{2\sigma^4}{n-p} \end{aligned}

証明

準備： 射影行列 \boldsymbol{M}_X の性質

行列 \boldsymbol{M}_X = \boldsymbol{I}_n - \boldsymbol{X} (\boldsymbol{X}^\top \boldsymbol{X})^{-1} \boldsymbol{X}^\top について、以下が成り立つ：

\begin{aligned} \boldsymbol{M}_X \boldsymbol{Y} = \hat{\boldsymbol{\varepsilon}},\\ \boldsymbol{M}_X \boldsymbol{X} = \boldsymbol{0}. \end{aligned}

\boldsymbol{M}_X は射影行列と呼ばれ、以下のような性質を持つ。

• 射影行列 \boldsymbol{M}_X は対称行列である。

\boldsymbol{M}_X^\top = \boldsymbol{M}_X

証明

\begin{aligned} \boldsymbol{M}_X^\top &= \left\{\boldsymbol{I}_n - \boldsymbol{X} (\boldsymbol{X}^\top \boldsymbol{X})^{-1} \boldsymbol{X}^\top \right\}^\top\\ &= \boldsymbol{I}_n - \boldsymbol{X} (\boldsymbol{X}^\top \boldsymbol{X})^{-1} \boldsymbol{X}^\top\\ &= \boldsymbol{M}_X \end{aligned}

• 射影行列 \boldsymbol{M}_X は冪等である。

\boldsymbol{M}_X^2 = \boldsymbol{M}_X

証明

\begin{aligned} \boldsymbol{M}_X^2 &= \left\{ \boldsymbol{I}_n - \boldsymbol{X} (\boldsymbol{X}^\top \boldsymbol{X})^{-1} \boldsymbol{X}^\top \right\} \left\{\boldsymbol{I}_n - \boldsymbol{X} (\boldsymbol{X}^\top \boldsymbol{X})^{-1} \boldsymbol{X}^\top \right\}\\ &= \boldsymbol{I}_n - 2\boldsymbol{X} (\boldsymbol{X}^\top \boldsymbol{X})^{-1} \boldsymbol{X}^\top \\ &\qquad + \boldsymbol{X} (\boldsymbol{X}^\top \boldsymbol{X})^{-1} \boldsymbol{X}^\top\boldsymbol{X} (\boldsymbol{X}^\top \boldsymbol{X})^{-1} \boldsymbol{X}^\top\\ &= \boldsymbol{I}_n - \boldsymbol{X} (\boldsymbol{X}^\top \boldsymbol{X})^{-1} \boldsymbol{X}^\top\\ &= \boldsymbol{M}_X \end{aligned}

• 射影行列 \boldsymbol{M}_X の階数は {\rm rank}\, \boldsymbol{M}_X = n-p である。

証明

\boldsymbol{X},\,\boldsymbol{M}_X の列ベクトルが張る空間をそれぞれ \mathcal{C}(\boldsymbol{X}),\, \mathcal{C}(\boldsymbol{M}_X) と表す。

{\rm rank} の定義より

{\rm rank}\, \boldsymbol{M}_X = {\rm dim}\, \mathcal{C}(\boldsymbol{M}_X)

と表せる。ここで、\mathcal{C}(\boldsymbol{M}_X) の直交補空間を \mathcal{C}(\boldsymbol{M}_X)^{\perp} とすると、

{\rm dim}\, \mathcal{C}(\boldsymbol{M}_X) + {\rm dim}\, \mathcal{C}(\boldsymbol{M}_X)^{\perp} = n

であり、さらに \boldsymbol{M}_X の核を {\rm Ker}\,\boldsymbol{M}_X = \left\{\boldsymbol{v}; \boldsymbol{M}_X \boldsymbol{v} = \boldsymbol{0} \right\} と表すと

\mathcal{C}(\boldsymbol{M}_X)^{\perp} = {\rm Ker}\,\boldsymbol{M}_X

が成り立つことから、^[1]

{\rm rank}\, \boldsymbol{M}_X = n - {\rm dim}\left({\rm Ker}\, \boldsymbol{M}_X\right)

のように表すことができる。

ここで、

{\rm Ker}\, \boldsymbol{M}_X = \mathcal{C}(\boldsymbol{X})

が成り立つ（証明は後述）。
これを用いると

{\rm dim}\left( {\rm Ker}\, \boldsymbol{M}_X \right)= {\rm dim}\, \mathcal{C}(\boldsymbol{X}) = {\rm rank}\, \boldsymbol{X} = p

が成り立つことから、

\begin{aligned} {\rm rank}\, \boldsymbol{M}_X &= n - {\rm dim}\left({\rm Ker}\, \boldsymbol{M}_X\right)\\ &= n-p \end{aligned}

が示された。

---

以降、{\rm Ker}\, \boldsymbol{M}_X = \mathcal{C}(\boldsymbol{X}) を示す。

\boldsymbol{M}_X \boldsymbol{X} = \boldsymbol{0} より、\mathcal{C}(\boldsymbol{X}) \subseteq {\rm Ker}\, \boldsymbol{M}_X は簡単に示すことができる。
なぜなら、 \boldsymbol{v} \in \mathcal{C}(\boldsymbol{X}) \Rightarrow \boldsymbol{M}_X\boldsymbol{v} = \boldsymbol{0} より \boldsymbol{v} \in {\rm Ker}\, \boldsymbol{M}_X となるからである。

次に {\rm Ker}\, \boldsymbol{M}_X \subseteq \mathcal{C}(\boldsymbol{X}) を示す。
そのためには、\boldsymbol{v} \notin \mathcal{C}(\boldsymbol{X}) \Rightarrow \boldsymbol{v} \notin {\rm Ker}\, \boldsymbol{M}_X が言えればよい。

\boldsymbol{v} \notin \mathcal{C}(\boldsymbol{X}) より、\boldsymbol{v} = \boldsymbol{v}_x + \boldsymbol{v}_{\bar{x}} \, \left( \boldsymbol{v}_x \in \mathcal{C}(\boldsymbol{X}), \, \boldsymbol{v}_{\bar{x}} \in \mathcal{C}(\boldsymbol{X})^\perp, \, \boldsymbol{v}_{\bar{x}} \ne \boldsymbol{0} \right) と表すことができる。
ここで、\boldsymbol{X}^\top \boldsymbol{v}_{\bar{x}} = \boldsymbol{0} より、

\begin{aligned} \boldsymbol{M}_X \boldsymbol{v}_{\bar{x}} &= \left\{\boldsymbol{I}_n - \boldsymbol{X} (\boldsymbol{X}^\top \boldsymbol{X})^{-1} \boldsymbol{X}^\top \right\}\boldsymbol{v}_{\bar{x}}\\ &=\boldsymbol{v}_{\bar{x}} \end{aligned}

が成り立つ。
従って、

\begin{aligned} \boldsymbol{M}_X \boldsymbol{v} &= \boldsymbol{M}_X \boldsymbol{v}_x + \boldsymbol{M}_X \boldsymbol{v}_{\bar{x}}\\ &=\boldsymbol{0} + \boldsymbol{v}_{\bar{x}}\\ &\ne \boldsymbol{0} \end{aligned}

となるため、\boldsymbol{v} \notin {\rm Ker}\, \boldsymbol{M}_X が示された。

以上より \mathcal{C}(\boldsymbol{X}) \subseteq {\rm Ker}\, \boldsymbol{M}_X かつ {\rm Ker}\, \boldsymbol{M}_X \subseteq \mathcal{C}(\boldsymbol{X}) より {\rm Ker}\, \boldsymbol{M}_X = \mathcal{C}(\boldsymbol{X}) が示された。

• 射影行列 \boldsymbol{M}_X は正規直交行列 \boldsymbol{P} を用いて固有値分解でき、さらにその固有値のうち n-p 個は1、残りの p 個は0となる。

\boldsymbol{M}_X = \boldsymbol{P}\boldsymbol{\Lambda}\boldsymbol{P}^\top

\boldsymbol{P}^\top \boldsymbol{P} = \boldsymbol{I}_n, \quad \boldsymbol{\Lambda} = {\rm diag}(\underbrace{1,1,...,1}_{n-p},\underbrace{0,0,...,0}_p)

証明

\boldsymbol{p}_i を \boldsymbol{M}_X の固有ベクトルとし、その固有値を \lambda_i とする。

このとき、

\boldsymbol{M}_X \boldsymbol{p}_i = \lambda_i \boldsymbol{p}_i

であり、さらに \boldsymbol{M}_X^2 = \boldsymbol{M}_X から

\begin{aligned} \boldsymbol{M}_X^2 \boldsymbol{p}_i &= \boldsymbol{M}_X \boldsymbol{p}_i = \lambda_i \boldsymbol{p}_i\\ &= \boldsymbol{M}_X \left( \boldsymbol{M}_X \boldsymbol{p}_i \right) = \boldsymbol{M}_X \left( \lambda_i \boldsymbol{p}_i\right) = \lambda_i^2 \boldsymbol{p}_i \end{aligned}

が成り立つため、固有値 \lambda_i は 0 or 1 であることがわかる。
さらに前述の {\rm rank}\, \boldsymbol{M}_X = n-p であることから、固有値のうち n-p 個が1で残りの p 個が0であることがわかる。

• 射影行列 \boldsymbol{M}_X の trace は {\rm tr}\, \boldsymbol{M}_X = n-p となる。

証明

先ほどのように \boldsymbol{M}_X = \boldsymbol{P} \boldsymbol{\Lambda} \boldsymbol{P}^\top と固有値分解すると、

\begin{aligned} {\rm tr}\, \boldsymbol{M}_X &= {\rm tr}\, \left( \boldsymbol{P} \boldsymbol{\Lambda} \boldsymbol{P}^\top \right)\\ &= {\rm tr}\, \left( \boldsymbol{\Lambda} \boldsymbol{P}^\top \boldsymbol{P} \right)\\ &= {\rm tr}\, \boldsymbol{\Lambda} \end{aligned}

と表せる。
対角行列 \Lambda の対角成分のうち n-p 個は 1 でそれ以外は 0 なので、{\rm tr}\, \boldsymbol{\Lambda} = n-p となる。
このことから、{\rm tr}\, \boldsymbol{M}_X = n-p が示された。

s^2 が誤差項分散の不偏推定量になることの証明

E[\hat{\boldsymbol{\varepsilon}}^\top\hat{\boldsymbol{\varepsilon}}] = (n-p)\sigma^2 を示す。

\begin{aligned} \hat{\boldsymbol{\varepsilon}} &= \boldsymbol{M}_X \boldsymbol{Y}\\ &= \boldsymbol{M}_X (\boldsymbol{X}\boldsymbol{\beta} + \boldsymbol{\varepsilon})\\ &= \boldsymbol{M}\boldsymbol{\varepsilon} \end{aligned}

であるので、

\begin{aligned} \hat{\boldsymbol{\varepsilon}}^\top\hat{\boldsymbol{\varepsilon}} &= (\boldsymbol{M}_X \boldsymbol{\varepsilon})^\top (\boldsymbol{M}_X \boldsymbol{\varepsilon})\\ &= \boldsymbol{\varepsilon}^\top \boldsymbol{M}_X \boldsymbol{\varepsilon} \end{aligned}

が成り立つ。
ここで、

\begin{aligned} \hat{\boldsymbol{\varepsilon}}^\top\hat{\boldsymbol{\varepsilon}} &={\rm tr}\,\left( \hat{\boldsymbol{\varepsilon}}^\top\hat{\boldsymbol{\varepsilon}} \right)\\ &= {\rm tr}\,\left( \boldsymbol{\varepsilon}^\top \boldsymbol{M}_X \boldsymbol{\varepsilon} \right)\\ &= {\rm tr}\,\left( \boldsymbol{M}_X \boldsymbol{\varepsilon} \boldsymbol{\varepsilon}^\top\right)\\ \end{aligned}

が成立し、さらに \boldsymbol{\varepsilon} \boldsymbol{\varepsilon}^\top = \sigma^2 \boldsymbol{I}_n また {\rm tr}\, \boldsymbol{M}_X = {\rm rank}\, \boldsymbol{M}_X = n-p であることから、

\begin{aligned} E\left[\hat{\boldsymbol{\varepsilon}}^\top\hat{\boldsymbol{\varepsilon}}\right] &=E\left[{\rm tr}\,\left( \boldsymbol{M}_X \boldsymbol{\varepsilon} \boldsymbol{\varepsilon}^\top\right)\right]\\ &= {\rm tr}\, \left( \boldsymbol{M}_X E[\boldsymbol{\varepsilon} \boldsymbol{\varepsilon}^\top ] \right)\\ &= \sigma^2 {\rm tr}\, \boldsymbol{M}_X\\ &= (n-p) \sigma^2 \end{aligned}

以上から、E[\hat{\boldsymbol{\varepsilon}}^\top\hat{\boldsymbol{\varepsilon}}] = (n-p)\sigma^2 が成り立つことが示された。

従って、

E[s^2] = E\left[ \frac{\hat{\boldsymbol{\varepsilon}}^\top\hat{\boldsymbol{\varepsilon}}}{n-p} \right] = \sigma^2

となる。

標準化した誤差項二乗和がカイ二乗分布に従うことの証明

次に、誤差項分散 \boldsymbol{\varepsilon} が独立同時な正規分布に従うとき \frac{\hat{\boldsymbol{\varepsilon}}^\top\hat{\boldsymbol{\varepsilon}}}{\sigma^2} \sim \chi^2(n-p) が成り立つことを示す。

\boldsymbol{M}_X は冪等行列なので、その固有値は、n-p 個が1で p 個が0である。
そこで、\boldsymbol{M}_X を固有値分解して以下のように表すことができる：

\boldsymbol{M}_X = \boldsymbol{P} \boldsymbol{\Lambda} \boldsymbol{P}^\top

ただし、\boldsymbol{P} は正規直交行列（\boldsymbol{P}^\top \boldsymbol{P} = \boldsymbol{I}_n）、\boldsymbol{\Lambda} は対角行列でその対角成分 \{\lambda_i\}_{i=1}^n は最初の n-p 個が1で残りの p 個が0である。

ここで、\hat{\boldsymbol{\varepsilon}}^\top \hat{\boldsymbol{\varepsilon}} は以下のように表すことができる：

\begin{aligned} \hat{\boldsymbol{\varepsilon}}^\top \hat{\boldsymbol{\varepsilon}} &= \boldsymbol{\varepsilon}^\top \boldsymbol{M}_X \boldsymbol{\varepsilon}\\ &= \boldsymbol{\varepsilon}^\top \boldsymbol{P} \boldsymbol{\Lambda} \boldsymbol{P}^\top \boldsymbol{\varepsilon}\\ &= \boldsymbol{e}^\top \boldsymbol{\Lambda} \boldsymbol{e}\\ &= \lambda_1 e_1^2 + \lambda_2 e_2^2 + ... + \lambda_{n-p}e_{n-p}^2 + \lambda_{n-p+1} e_{n-p+1}^2 + ... + \lambda_n e_n^2\\ &= e_1^2 + e_2^2 + ... + e_{n-p}^2 \end{aligned}

途中、\boldsymbol{e} = \boldsymbol{P}^\top \boldsymbol{\varepsilon} と置いた。
E[\boldsymbol{e}] = \boldsymbol{0}, V[\boldsymbol{e}] = \sigma^2 \boldsymbol{P}^\top \boldsymbol{I}_n \boldsymbol{P} = \sigma^2 \boldsymbol{I}_n であることから、 \boldsymbol{e} \sim N(\boldsymbol{0}, \sigma^2 \boldsymbol{I}_n) である。
したがって、e_i / \sigma \sim N(0,1) かつ \{e_i\}_{i=1}^n は独立同時。

このことから、

e_1^2 / \sigma^2 + e_2^2 /\sigma^2 + ... + e_{n-p}^2/\sigma^2 \sim \chi^2(n-p)

であり、以上から

\frac{\hat{\boldsymbol{\varepsilon}}^\top \hat{\boldsymbol{\varepsilon}}}{\sigma^2} \sim \chi^2(n-p)

が成立することがわかる。

参考文献

• 佐和隆光、回帰分析(新装版) (統計ライブラリー、朝倉書店、2020)

https://amzn.asia/d/2vsEmZ6

• J.D.Hamilton, Time Series Analysis (Princeton Univ Pr, 1994)

https://amzn.asia/d/dHAwb2G

脚注

1. https://zenn.dev/link/comments/116bc995370fff の補足2参照 ↩︎