🔔 本の下書きプレビューをシェアできるようになりました
Zenn
Open6

線形回帰の復習

統計学
計量経済学
回帰分析
線形回帰
畳屋民也畳屋民也

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

設定

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

yi=β0+β1x1i+β2x2i+...βpxp−1,i+εi=xi⊤β+εi. \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}

ここで、

xi=(1x1ix2i⋮xp−1,i),β=(β0β1β2⋮βp−1) \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}

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

E[εi]=0,E[εiεj]={σ2i=j0i≠j. 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. .

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

NN個の観測値{(x1,y1),(x2,y2),...,(xN,yN)}\{ (\boldsymbol{x}_1, y_1), (\boldsymbol{x}_2, y_2), ..., (\boldsymbol{x}_N, y_N) \} が与えられているとき、

Y=(y1y2⋮yN),ε=(ε1ε2⋮ε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},
X=(x1⊤x2⊤⋮xN⊤)=(1x11x21...xp−1,11x12x22...xp−1,2⋮⋮⋮1x1Nx2N...xp−1,N) \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}

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

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

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

yi=β0+β1xi+εi y_i = \beta_0 + \beta_1 x_i + \varepsilon_i

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

X=(1x11x2⋮1xN),β=(β0β1) \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}

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

(y1y2⋮yN)=(1x11x2⋮1xN)(β0β1)+(ε1ε2⋮εN)=(β0+β1x1+ε1β0+β1x2+ε2⋮β0+β1xN+εN) \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}
畳屋民也畳屋民也

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

回帰式

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

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

β^=(X⊤X)−1X⊤Y=(∑i=1Nxixi⊤)−1(∑i=1Nxiyi)(1) \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}

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

E[β^]=β.(2) E[\hat{\boldsymbol{\beta}}] = \boldsymbol{\beta}. \tag{2}

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

V[β^]=σ2(X⊤X)−1=σ2(∑i=1Nxixi⊤)−1.(3) \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}

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

β^∼N(β,σ2(X⊤X)−1) \hat{\boldsymbol{\beta}} \sim N(\boldsymbol{\beta}, \sigma^2 (\boldsymbol{X}^\top \boldsymbol{X})^{-1})

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

yi=β0+β1xi+εi y_i = \beta_0 + \beta_1 x_i + \varepsilon_i

の場合を考える。

係数推定量

このとき、

X=(1x11x2⋮1xN),β=(β0β1) \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}

であるので、

(X⊤X)−1=(∑i1∑ixi∑ixi∑ixi2)−1=1N∑ixi2−(∑ixi)2(∑ixi2−∑ixi−∑ixiN)=1∑i(xi−xˉ)2(1N∑ixi2−xˉ−xˉ1) \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}

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

β^=1∑i(xi−xˉ)2(1N∑ixi2−xˉ−xˉ1)(∑iyi∑ixiyi) \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}

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

これを整理すると、

β^0=yˉ∑i=1Nxi2−xˉ∑i=1Nxiyi∑i=1N(xi−xˉ)2=yˉ−β^1xˉβ^1=∑i=1N(yi−yˉ)(xi−xˉ)∑i=1N(xi−xˉ)2 \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}

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

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

前述のように

(X⊤X)−1=1∑i(xi−xˉ)2(1N∑ixi2−xˉ−xˉ1) (\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[β^]=σ2(X⊤X)−1V[\hat{\boldsymbol{\beta}}] = \sigma^2 (\boldsymbol{X}^\top \boldsymbol{X})^{-1} に代入すると、

V[β^0]=σ21N∑i=1Nxi2∑i=1N(xi−xˉ)2V[β^1]=σ2∑i=1N(xi−xˉ)2Cov[β^0,β^1]=−σ2xˉ∑i=1N(xi−xˉ)2 \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}

係数推定量の分布

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

(c^ϕ^)∼N((cϕ),σ2∑i(xi−xˉ)2(1N∑ixi2−xˉ−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) の導出

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

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

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

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

∂RSS(β)∂β=∂(Y⊤Y−Y⊤Xβ−β⊤X⊤Y+β⊤X⊤Xβ)∂β=∂(Y⊤Y)∂β−2∂(Y⊤Xβ)∂β+∂(β⊤X⊤Xβ)∂β=0−2(Y⊤X)⊤+{X⊤X+(X⊤X)⊤}β=−2X⊤Y+2X⊤Xβ \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行目では、Y⊤Xβ\boldsymbol{Y}^\top \boldsymbol{X}\boldsymbol{\beta} がスカラーであることからβ⊤X⊤Y=(Y⊤Xβ)⊤=Y⊤Xβ\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行目では、以下のベクトルによる微分の関係式を適用した：

∂(Ax)∂x=A⊤,∂(x⊤Ax)∂x=(A+A⊤)x. \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(β)RSS(\boldsymbol{\beta}) が最小になるとき上記の微分の値は0\boldsymbol{0} になるので、

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

より、

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

を解いて、

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

(2) の証明

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

β^=(X⊤X)−1X⊤Y=(X⊤X)−1X⊤(Xβ+ε)=(X⊤X)−1X⊤Xβ+(X⊤X)−1X⊤ε=β+(X⊤X)−1X⊤ε(∵(X⊤X)−1(X⊤X)=1) \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}

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

E[β^]=β+(X⊤X)−1X⊤E[ε]=β+0=β \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[β^]=σ2(X⊤X)−1V[\hat{\boldsymbol{\beta}}] = \sigma^2 (\boldsymbol{X}^\top \boldsymbol{X})^{-1} を示す。

V[β^]=E[(β^−β)(β^−β)⊤]=E[{(X⊤X)−1X⊤ε}{(X⊤X)−1X⊤ε}⊤]=(X⊤X)−1X⊤E[εε⊤](X⊤)⊤{(X⊤X)−1}⊤=(X⊤X)−1X⊤σ21X(X⊤X)−1=σ2(X⊤X)−1(X⊤X)(X⊤X)−1=σ2(X⊤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}

補足

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

xi=(1x1ix2i⋮xpi),X=(x1⊤x2⊤⋮xN⊤)=(1x11x21...xp11x12x22...xp2⋮⋮⋮1x1Nx2N...xpN) \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}

であるので、

X⊤X=∑i=1Nxixi⊤=(∑i1∑ix1i∑ix2i...∑ixpi∑ix1i∑ix1i2∑ix1ix2i...∑ix1ixpi∑ix2i∑ix2ix1i∑ix2i2...∑ix2ixpi⋮⋮⋮∑ixpi∑ixpix1i∑ixpix2i...∑ixpi2) \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 回帰とその推定量の性質

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

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

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

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

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

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

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

結果

サンプルサイズを nn、説明変数の数を pp とする。
このとき、X∈Rn×p, Y∈Rn×1, β∈Rp×1\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} の推定量は、以下である：

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

IpI_p は p×pp \times p の単位行列である。

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

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

Ridge 回帰の概要

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

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

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

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

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

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

∂L(β)∂β=−2X⊤(Y−Xβ)+2αβ=−2X⊤Y+2(X⊤X+αIp)β. \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(β)L(\boldsymbol{\beta}) は最小になるβ=β^α\boldsymbol{\beta} = \hat{\boldsymbol{\beta}}_\alpha では、上記の微分の値が 0\boldsymbol{0} になる。

従ってこのとき、

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

が成り立つ。

ここで、行列 (X⊤X+αIp)(\boldsymbol{X}^\top \boldsymbol{X} + \alpha I_p) は α>0\alpha > 0 で常に非特異で逆行列を持つ[1]ので、先ほどの式の両辺に左から (X⊤X+αIp)−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[β^α]≠βE[\hat{\boldsymbol{\beta}}_\alpha] \ne \boldsymbol{\beta}）ことを示す。

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

β^α=(X⊤X+αIp)−1X⊤(Xβ+ε)=(X⊤X+αIp)−1X⊤Xβ+(X⊤X+αIp)−1X⊤ε(3) \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}

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

E[β^α]=(X⊤X+αIp)−1X⊤Xβ=(X⊤X+αIp)−1(X⊤X+αIp)β−(X⊤X+αIp)−1αIpβ=β−α(X⊤X+αIp)−1β. \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} と一致せず、−α(X⊤X+αIp)−1β- \alpha (\boldsymbol{X}^\top \boldsymbol{X} + \alpha I_p)^{-1} \boldsymbol{\beta} だけバイアスが入る。

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

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

式(3) より、

V[β^α]=E[{(X⊤X+αIp)−1X⊤ε}{(X⊤X+αIp)−1X⊤ε}⊤]=(X⊤X+αIp)−1X⊤E[εε⊤]X{(X⊤X+αIp)−1}⊤=σ2(X⊤X+αIp)−1X⊤X(X⊤X+αIp)−1 \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[εε⊤]=σ2In E[\boldsymbol{\varepsilon} \boldsymbol{\varepsilon}^\top] = \sigma^2 I_n
{(X⊤X+αIp)−1}⊤=(X⊤X+αIp)−1 \{ (\boldsymbol{X}^\top \boldsymbol{X} + \alpha I_p)^{-1} \}^\top = (\boldsymbol{X}^\top \boldsymbol{X} + \alpha I_p)^{-1}

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

脚注

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

Hidden comment
Hidden comment
畳屋民也畳屋民也

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

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

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

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

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

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

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

ε^=Y−Xβ^={In−X(X⊤X)−1X⊤}Y. \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}

以降、 MX=In−X(X⊤X)−1X⊤\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}} について、

s2=ε^⊤ε^n−p s^2 = \frac{\hat{\boldsymbol{\varepsilon}}^\top\hat{\boldsymbol{\varepsilon}}}{n-p}

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

E[s2]=σ2 E[s^2] = \sigma^2

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

ε^⊤ε^σ2∼χ2(n−p) \frac{\hat{\boldsymbol{\varepsilon}}^\top\hat{\boldsymbol{\varepsilon}}}{\sigma^2} \sim \chi^2(n-p)

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

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

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

証明

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

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

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

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

  • 射影行列 MX\boldsymbol{M}_X は対称行列である。
MX⊤=MX \boldsymbol{M}_X^\top = \boldsymbol{M}_X
証明
MX⊤={In−X(X⊤X)−1X⊤}⊤=In−X(X⊤X)−1X⊤=MX \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}
  • 射影行列 MX\boldsymbol{M}_X は冪等である。
MX2=MX \boldsymbol{M}_X^2 = \boldsymbol{M}_X
証明
MX2={In−X(X⊤X)−1X⊤}{In−X(X⊤X)−1X⊤}=In−2X(X⊤X)−1X⊤+X(X⊤X)−1X⊤X(X⊤X)−1X⊤=In−X(X⊤X)−1X⊤=MX \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}
  • 射影行列 MX\boldsymbol{M}_X の階数は rank MX=n−p{\rm rank}\, \boldsymbol{M}_X = n-p である。
証明

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

rank{\rm rank} の定義より

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

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

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

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

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

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

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

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

ここで、

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

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

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

が成り立つことから、

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

が示された。

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

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

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

v∉C(X)\boldsymbol{v} \notin \mathcal{C}(\boldsymbol{X}) より、v=vx+vxˉ (vx∈C(X), vxˉ∈C(X)⊥, vxˉ≠0)\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) と表すことができる。
ここで、X⊤vxˉ=0\boldsymbol{X}^\top \boldsymbol{v}_{\bar{x}} = \boldsymbol{0} より、

MXvxˉ={In−X(X⊤X)−1X⊤}vxˉ=vxˉ \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}

が成り立つ。
従って、

MXv=MXvx+MXvxˉ=0+vxˉ≠0 \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}

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

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

  • 射影行列 MX\boldsymbol{M}_X は正規直交行列 P\boldsymbol{P} を用いて固有値分解でき、さらにその固有値のうち n−pn-p 個は1、残りの pp 個は0となる。
MX=PΛP⊤ \boldsymbol{M}_X = \boldsymbol{P}\boldsymbol{\Lambda}\boldsymbol{P}^\top
P⊤P=In,Λ=diag(1,1,...,1⏟n−p,0,0,...,0⏟p) \boldsymbol{P}^\top \boldsymbol{P} = \boldsymbol{I}_n, \quad \boldsymbol{\Lambda} = {\rm diag}(\underbrace{1,1,...,1}_{n-p},\underbrace{0,0,...,0}_p)
証明

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

このとき、

MXpi=λipi \boldsymbol{M}_X \boldsymbol{p}_i = \lambda_i \boldsymbol{p}_i

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

MX2pi=MXpi=λipi=MX(MXpi)=MX(λipi)=λi2pi \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}

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

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

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

tr MX=tr (PΛP⊤)=tr (ΛP⊤P)=tr Λ \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−pn-p 個は 1 でそれ以外は 0 なので、tr Λ=n−p{\rm tr}\, \boldsymbol{\Lambda} = n-p となる。
このことから、tr MX=n−p{\rm tr}\, \boldsymbol{M}_X = n-p が示された。

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

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

ε^=MXY=MX(Xβ+ε)=Mε \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}

であるので、

ε^⊤ε^=(MXε)⊤(MXε)=ε⊤MXε \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}

が成り立つ。
ここで、

ε^⊤ε^=tr (ε^⊤ε^)=tr (ε⊤MXε)=tr (MXεε⊤) \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}

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

E[ε^⊤ε^]=E[tr (MXεε⊤)]=tr (MXE[εε⊤])=σ2tr MX=(n−p)σ2 \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[ε^⊤ε^]=(n−p)σ2E[\hat{\boldsymbol{\varepsilon}}^\top\hat{\boldsymbol{\varepsilon}}] = (n-p)\sigma^2 が成り立つことが示された。

従って、

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

となる。

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

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

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

MX=PΛP⊤ \boldsymbol{M}_X = \boldsymbol{P} \boldsymbol{\Lambda} \boldsymbol{P}^\top

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

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

ε^⊤ε^=ε⊤MXε=ε⊤PΛP⊤ε=e⊤Λe=λ1e12+λ2e22+...+λn−pen−p2+λn−p+1en−p+12+...+λnen2=e12+e22+...+en−p2 \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}

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

このことから、

e12/σ2+e22/σ2+...+en−p2/σ2∼χ2(n−p) e_1^2 / \sigma^2 + e_2^2 /\sigma^2 + ... + e_{n-p}^2/\sigma^2 \sim \chi^2(n-p)

であり、以上から

ε^⊤ε^σ2∼χ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参照 ↩︎

ログインするとコメントできます