線形回帰分析の残差の不偏分散 \(s^2\) を求めます。
参考にしました資料は https://user.keio.ac.jp/~nagakura/zemi/K_regression.pdf です。
以下の多変量回帰分析(切片項あり)を考えます。 \[\mathbf{y}=\mathbf{X}\boldsymbol{\beta}+\boldsymbol{\epsilon}\] ここで、 \[\mathbf{y}=\begin{bmatrix}y_1\\y_2\\\vdots\\y_n\end{bmatrix}\quad\mathbf{X}=\begin{bmatrix}1&x_{21} & \cdots& x_{k1} \\1 & x_{22}&\cdots& x_{k2}\\\vdots&\vdots&\ddots&\vdots\\1&x_{2n}&\cdots&x_{kn}\end{bmatrix}\quad\boldsymbol{\epsilon}=\begin{bmatrix}\epsilon_1\\\epsilon_2\\\vdots\\\epsilon_n\end{bmatrix}\quad\boldsymbol{\beta}=\begin{bmatrix}\beta_1\\\beta_2\\\vdots\\\beta_k\end{bmatrix}\] \(\boldsymbol{\beta}\) のOLS推定量を \(\mathbf{b}=\left(b_1,b_2\cdots,b_k\right)^{T}\) とすると、残差平方和は以下の通りとなります。 \[\displaystyle\sum_{i=1}^n\left(y_i-b_1-b_2x_{2i}-\cdots-b_kx_{ki}\right)^2=\left(\mathbf{y}-\mathbf{Xb}\right)^{T}\left(\mathbf{y}-\mathbf{Xb}\right)\] 残差平方和 \(\left(\mathbf{y}-\mathbf{Xb}\right)^{T}\left(\mathbf{y}-\mathbf{Xb}\right)\) を \(\mathbf{b}\) に関して偏微分した結果をゼロとおき、 \[-2\mathbf{X}^{T}\mathbf{y}+-2\mathbf{X}^{T}\mathbf{Xb}=0\] \(\mathbf{b}\) について解いた、 \[\mathbf{b}=\left(\mathbf{X}^{T}\mathbf{X}\right)^{-1}\mathbf{X}^{T}\mathbf{y}\] が、\(\boldsymbol{\beta}\) のOLS推定量となります(以降、\(\boldsymbol{\hat{\beta}}\) と表記します)。
続いて誤差項 \(\epsilon\) の分散 \(\sigma^2\) の不偏推定量 \(s^2\) を求めます。
\(\epsilon\) の推定量を \(\mathbf{e}=\left(e_1,e_2,\cdots,e_n\right)^T\) とすると、 \[\begin{eqnarray}\mathbf{e}&=&\mathbf{y}-\mathbf{X}\boldsymbol{\hat{\beta}}\\&=&\mathbf{y}-\mathbf{X}\left(\mathbf{X}^{T}\mathbf{X}\right)^{-1}\mathbf{X}^{T}\mathbf{y}\\&=&\left(\mathbf{I}-\mathbf{X}\left(\mathbf{X}^{T}\mathbf{X}\right)^{-1}\mathbf{X}^{T}\right)\mathbf{y}\end{eqnarray}\] ここで、\(\mathbf{M}=\mathbf{I}-\mathbf{X}\left(\mathbf{X}^{T}\mathbf{X}\right)^{-1}\mathbf{X}^{T}\) とすると、 \[\mathbf{e}=\mathbf{My}\] さらに、 \[\begin{eqnarray}\mathbf{MX}&=&\left(\mathbf{I}-\mathbf{X}\left(\mathbf{X}^{T}\mathbf{X}\right)^{-1}\mathbf{X}^{T}\right)\mathbf{X}\\&=&\left(\mathbf{X}-\mathbf{X}\left(\mathbf{X}^{T}\mathbf{X}\right)^{-1}\mathbf{X}^{T}\mathbf{X}\right)\\&=&\mathbf{X}-\mathbf{XI}\\&=&\mathbf{X}-\mathbf{X}\\&=&0\end{eqnarray}\] であるので、 \[\begin{eqnarray}\mathbf{e}=\mathbf{My}=\mathbf{M}\left(\mathbf{X}\boldsymbol{\beta}+\boldsymbol{\epsilon}\right)=\mathbf{MX}\boldsymbol{\beta}+\mathbf{M}\boldsymbol{\epsilon}=\mathbf{M}\boldsymbol{\epsilon}\end{eqnarray}\] また、\(\mathbf{MM}\) は、 \[\begin{eqnarray} \mathbf{M}\mathbf{M} &=&\left(\mathbf{I}-\mathbf{X}\left(\mathbf{X}^{T}\mathbf{X}\right)^{-1}\mathbf{X}^{T}\right)\left(\mathbf{I}-\mathbf{X}\left(\mathbf{X}^{T}\mathbf{X}\right)^{-1}\mathbf{X}^{T}\right)\\ &=&\mathbf{I}^2-\mathbf{I}\left(\mathbf{X}\left(\mathbf{X}^{T}\mathbf{X}\right)^{-1}\mathbf{X}^{T}\right)-\left(\mathbf{X}\left(\mathbf{X}^{T}\mathbf{X}\right)^{-1}\mathbf{X}^{T}\right)\mathbf{I}+\left(\mathbf{X}\left(\mathbf{X}^{T}\mathbf{X}\right)^{-1}\mathbf{X}^{T}\right)^2\\ &=&\mathbf{I}-2\left(\mathbf{X}\left(\mathbf{X}^{T}\mathbf{X}\right)^{-1}\mathbf{X}^{T}\right)+\left(\mathbf{X}\left(\mathbf{X}^{T}\mathbf{X}\right)^{-1}\mathbf{X}^{T}\right)\left(\mathbf{X}\left(\mathbf{X}^{T}\mathbf{X}\right)^{-1}\mathbf{X}^{T}\right)\\ &=&\mathbf{I}-2\left(\mathbf{X}\left(\mathbf{X}^{T}\mathbf{X}\right)^{-1}\mathbf{X}^{T}\right)+\mathbf{X}\left(\mathbf{X}^{T}\mathbf{X}\right)^{-1}\left(\mathbf{X}^{T}\mathbf{X}\right)\left(\mathbf{X}^{T}\mathbf{X}\right)^{-1}\mathbf{X}^{T}\\ &=&\mathbf{I}-2\left(\mathbf{X}\left(\mathbf{X}^{T}\mathbf{X}\right)^{-1}\mathbf{X}^{T}\right)+\mathbf{X}\left(\mathbf{X}^{T}\mathbf{X}\right)^{-1}\mathbf{I}\mathbf{X}^{T}\\ &=&\mathbf{I}-2\left(\mathbf{X}\left(\mathbf{X}^{T}\mathbf{X}\right)^{-1}\mathbf{X}^{T}\right)+\mathbf{X}\left(\mathbf{X}^{T}\mathbf{X}\right)^{-1}\mathbf{X}^{T}\\ &=&\mathbf{I}-\mathbf{X}\left(\mathbf{X}^{T}\mathbf{X}\right)^{-1}\mathbf{X}^{T}\\ &=&\mathbf{M} \end{eqnarray}\] となることから、\(\mathbf{M}\) は冪等行列(\(\mathbf{M}\mathbf{M}=\mathbf{M}\))であり、さらに \(\mathbf{M}\)の転置行列は、 \[\begin{eqnarray} \mathbf{M}^T&=&\left(\mathbf{I}-\mathbf{X}\left(\mathbf{X}^{T}\mathbf{X}\right)^{-1}\mathbf{X}^{T}\right)^T\\ &=&\mathbf{I}^T-\left(\mathbf{X}\left(\mathbf{X}^{T}\mathbf{X}\right)^{-1}\mathbf{X}^{T}\right)^T\\ &=&\mathbf{I}-\left(\mathbf{X}^T\right)^T\left(\left(\mathbf{X}^{T}\mathbf{X}\right)^{-1}\right)^T\mathbf{X}^T\\ &=&\mathbf{I}-\mathbf{X}\left(\left(\mathbf{X}^{T}\mathbf{X}\right)^T\right)^{-1}\mathbf{X}^T\\ &=&\mathbf{I}-\mathbf{X}\left(\mathbf{X}^{T}\left(\mathbf{X}^T\right)^T\right)^{-1}\mathbf{X}^T\\ &=&\mathbf{I}-\mathbf{X}\left(\mathbf{X}^{T}\mathbf{X}\right)^{-1}\mathbf{X}^T\\ &=&\mathbf{M} \end{eqnarray}\] となることから、\(\mathbf{M}\) は対称行列(\(\mathbf{M}^T=\mathbf{M}\))でもあるため、 \[ \begin{eqnarray} \displaystyle\sum_{i=1}^ne_i^2&=&\mathbf{e}^T\mathbf{e}=\left(\mathbf{M}\boldsymbol{\epsilon}\right)^T\left(\mathbf{M}\boldsymbol{\epsilon}\right)=\boldsymbol{\epsilon}^T\left(\mathbf{M}^T\mathbf{M}\right)\boldsymbol{\epsilon}=\boldsymbol{\epsilon}^T\left(\mathbf{M}\mathbf{M}\right)\boldsymbol{\epsilon}=\boldsymbol{\epsilon}^T\mathbf{M}\boldsymbol{\epsilon} \end{eqnarray} \] となります。
ここで、\(\mathrm{E}\left(\epsilon_i\right)=0\quad i=1,\cdots,n\) かつ \(\mathrm{Cov}\left(\epsilon_i,\epsilon_j\right)=0\quad i\neq j,i,j=1,\cdots,n\) と仮定しているため \[\mathrm{Cov}\left(\epsilon_i,\epsilon_j\right)=\mathrm{E}\left(\left(\epsilon_i-\mathrm{E}\left(\epsilon_i\right)\right)\left(\epsilon_j-\mathrm{E}\left(\epsilon_j\right)\right)\right)=\mathrm{E}\left(\epsilon_i\epsilon_j\right)=0\quad i\neq j\]となり、よって \[\displaystyle{\sum_{i=1}^n}\displaystyle{\sum_{j=1}^n}\mathrm{E}\left(\epsilon_i\epsilon_j\right)=\displaystyle\sum_{i=1}^n\mathrm{E}\left(\epsilon_{ii}\right)=\sigma^2\] であることから、その期待値は \[\mathrm{E}\left(\boldsymbol{\epsilon}^T\mathbf{M}\boldsymbol{\epsilon}\right)=\mathrm{E}\left(\displaystyle\sum_{i=1}^n\displaystyle\sum_{j=1}^nm_{ij}\epsilon_i\epsilon_j\right)=\displaystyle\sum_{i=1}^n\displaystyle\sum_{j=1}^nm_{ij}\mathrm{E}\left(\epsilon_i\epsilon_j\right)=\sigma^2\displaystyle\sum_{i=1}^nm_{ii}=\sigma^2\,\mathrm{tr}(\mathbf{M}) \] ここで、\(m_{ij}\) は \(\mathbf{M}\) の \((i,j)\) 成分であり、\(\mathrm{tr}\left(\mathbf{M}\right)\) は行列 \(\mathbf{M}\) の対角要素の和を表します。
前述の通り \(\mathbf{M}\) は対称冪等行列であるので、そのランクと対角要素の和は一致します(参照資料: https://ymurasawa.web.fc2.com/ue-sl11.pdf)。
よって、 \[\begin{eqnarray} \mathrm{tr}\left(\mathbf{M}\right)&=&\mathrm{tr}\left(\mathbf{I}-\mathbf{X}\left(\mathbf{X}^{T}\mathbf{X}\right)^{-1}\mathbf{X}^{T}\right)\\ &=&\mathrm{tr}\left(\mathbf{I}\right)-\mathrm{tr}\left(\mathbf{X}\left(\mathbf{X}^{T}\mathbf{X}\right)^{-1}\mathbf{X}^{T}\right)\\ &=&\mathrm{rk}\left(\mathbf{I}\right)-\mathrm{rk}\left(\mathbf{X}\left(\mathbf{X}^{T}\mathbf{X}\right)^{-1}\mathbf{X}^{T}\right)\\ &=&n-k\end{eqnarray}\] となり、 \[\mathrm{E}\left(\boldsymbol{\epsilon}^T\mathbf{M}\boldsymbol{\epsilon}\right)=\sigma^2(n-k)\] よって \(s^2\) の期待値は \[\begin{eqnarray}\mathrm{E}\left(s^2\right)=\sigma^2&=&\dfrac{1}{n-k}\mathrm{E}\left(\boldsymbol{\epsilon}^T\mathbf{M}\boldsymbol{\epsilon}\right)\\ &=&\dfrac{1}{n-k} \displaystyle\sum_{i=1}^ne_i^2\end{eqnarray}\]
以下のサンプルを最小二乗法により解き、不偏分散を求めます。
始めに偏回帰係数を求めます。
set.seed(20240628)
n <- 30
x1 <- runif(n = n)
x2 <- runif(n = n)
b <- 1
X <- as.matrix(cbind(b, x1, x2))
y <- 2 * x1 + 4 * x2 + b + rnorm(n = n)
beta <- solve(t(X) %*% X) %*% t(X) %*% y
beta [,1]
b 0.2530877
x1 2.7824476
x2 5.2335449続いて残差の標準誤差(\(\sqrt{\sigma^2}\))を求めます。
lm.residual <- y - X %*% beta
n <- nrow(X)
k <- ncol(X)
sqrt(sum(lm.residual^2) / (n - k))[1] 1.12464関数 lm {stats} の結果と比較します。
library(dplyr)
lm(y ~ x1 + x2) %>% summary()
Call:
lm(formula = y ~ x1 + x2)
Residuals:
Min 1Q Median 3Q Max
-2.97617 -0.55784 0.00417 0.68928 2.24403
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 0.2531 0.5613 0.451 0.65565
x1 2.7824 0.7758 3.587 0.00131 **
x2 5.2335 0.8193 6.388 7.66e-07 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 1.125 on 27 degrees of freedom
Multiple R-squared: 0.68, Adjusted R-squared: 0.6563
F-statistic: 28.69 on 2 and 27 DF, p-value: 2.086e-07残差の標準誤差( Residual standard error )が一致しています。
以上です。
