5.1 線形回帰モデル
5.1.1 線形単回帰モデル
目的変数$y$の期待値は説明変数$x_i$の関数$E[y]=f(x_i)$で表され、誤差項$ε_i$は互いに独立で、正規分布$N(0,σ^2)$に従うとすると次の式で表すことができる。$y_i$の期待値は$x_i$の関数$E[y_i|x_i]=f(x_i)$で表される。
y_i=f(x_i)+ε_i, \quad ε\sim N(0,σ^2) \quad (i=1,\dots, n)
単回帰モデルでは次の通りになる。
y_i=α+βx_i+ε_i, \quad ε\sim N(0,σ^2) \quad (i=1,\dots, n)
確率変数$x,y$とすると、説明変数$x$を与えた時の$y$の条件付き期待値は$E[y|x]=α+βx$、分散は$V[y|x]=σ^2$で、誤差項は互いに独立であり、$x$とは独立に正規分布に従うと仮定する。
5.1.2 回帰係数の区間推定
単回帰モデルでは、最小二乗法によって得られた$\hat{β}$は確率変数であり、$β$の不偏推定量である。
\hat{\beta}_1 = \frac{\sum_{i=1}^n (x_i - \bar{x})(y_i - \bar{y})}
{\sum_{i=1}^n (x_i - \bar{x})^2}
= \boxed{\frac{T_{xy}}{T_{xx}}}
このことを確かめるために$T_{xy}=\sum_{i=1}^{n}(x_i-\bar{x})y_i$および$\sum_{i=1}^{n}(x_i-\bar{x})=0$という性質を利用し、$y_i$に$α+βx_i+ε_i$を代入して整理すると、次の式が得られる。最後の式の第2項に確率変数$ε_i$があるため、$\hat{β}$が確率変数であることが分かる。
T_{xy}=\sum_{i=1}^{n}(x_i-\bar{x})(α+βx_i+ε_i)=β\sum_{i=1}^{n}(x_i-\bar{x})x_i+\sum_{i=1}^{n}(x_i-\bar{x})ε_i
また、$\sum_{i=1}^{n}(x_i-\bar{x})x_i=\sum_{i=1}^{n}[(x_i-\bar{x})^2-(x_i-\bar{x})\bar{x}]=T_{xx}+0=T_{xx}$であるから
=βT_{xx}+\sum_{i=1}^{n}(x_i-\bar{x})ε_i
なお、 $T_{xy} = βT_{xx}+\sum_{i=1}^{n}(x_i-\bar{x})ε_i$ が成り立つ証明は次の通り。
$T_{xy}$ の定義は次の通り：
T_{xy} = \sum_{i=1}^{n}(x_i - \bar{x})(y_i - \bar{y})= \sum_{i=1}^{n}(x_i - \bar{x})y_i - \sum_{i=1}^{n}(x_i - \bar{x})\bar{y}
ここで $\bar{y}$ は定数なので外へ出せるので
= \sum_{i=1}^{n}(x_i - \bar{x})y_i - \bar{y}\sum_{i=1}^{n}(x_i - \bar{x})
残差は平方和ではないのであれば$0$となる性質より、$\sum_{i=1}^{n}(x_i - \bar{x}) = 0$となるため
第2項は 0 となる。よって次の式が売られる。
T_{xy} = \sum_{i=1}^{n}(x_i - \bar{x})y_i
また、あるデータ$x_i$に対する残差平方和への寄与率$w_i=(x_i-\bar{x})/T_{xx}$として、下式を$T_{xx}$で割ると
T_{xy}=βT_{xx}+\sum_{i=1}^{n}(x_i-\bar{x})ε_i
\frac{T_{xy}}{T_{xx}}=\hat{β}=β+\sum_{i=1}^{n}\frac{(x_i-\bar{x})ε_i}{T_{xx}}=β+\sum_{i=1}^{n}w_iε_i
$ε_i$が正規分布に従い、なおかつ上式は一次結合の形であるので、$\hat{β}$も正規分布に従うことがわかる。
また、$\hat{β}$の期待値と分散は次の通り。$x$は確率変数ではなく固定値として扱うため、次の通りの計算となる。
E[\hat{β}]=β+\sum_{i=1}^{n}w_iE[ε_i]=β
V[\hat{β}]=E[(\hat{β}-β)^2]=E[(\sum_{i=1}^{n}w_iε_i)^2]=\sum_{i=1}^{n}w_i^2E[ε_i^2]=\sum_{i=1}^{n}w_i^2Var[ε_i]
=σ^2*\sum_{i=1}^{n}w_i^2=σ^2*\sum_{i=1}^{n}\frac{(x_i-\bar{x})^2}{T_{xx}^2}=σ^2*\frac{T_{xx}}{T_{xx}^2}=\frac{σ^2}{T_{xx}}
($E[(\sum_{i=1}^{n}w_iε_i)^2]$を展開すると、$2\sum_{1≦i<j≦n} w_iw_jε_iε_j$の交差項はでてくるが独立性と$E[ε_i]=0$より消える。)
よって、$\hat{β}$は次の通りの正規分布に従う。
\boxed{\hat{β} \sim N(β,\frac{σ^2}{T_{xx}})}
誤差項$ε_i$の分散$σ^2$を推定するには、理論的に$y$の予測値と実測値の残差$(\frac{\hat{y_i}-y_i}{n})$が推定量を構成すると考えられる。この場合、最小二乗法でパラメーターを2つを求めている理論構成のため自由度を2を消費しているため、不偏推定量は次の通り。
\boxed{\hat{σ}^2=\frac{\sum_{i=1}^n(\hat{y_i}-y_i)^2}{n-2}=\frac{T_{yy}}{n-2}}
さて、ようやく$\hat{β}$の信頼区間を導出できる。ここで、統計検定量$t$は次の通りで、$n-2$の自由度のt分布に従うことで95%信頼区間を算出できる。
\boxed{t=\frac{\hat{β}-β}{σ/\sqrt{T_{xx}} } \sim t_{(n-2)}}
β=\hat{β}±t_{0.025(n-2)}×σ/\sqrt{T_{xx}}
手計算の時には以下の式を使用すると良い(試験対策)。
\sum_{i=1}^n(\hat{y_i}-y_i)^2=T_{yy}-\hat{β}T_{xy}
- 以降は教科書の範囲外(切片の分散と回帰直線の予測値の信頼区間)
同様に切片の分散を導出できる。回帰式は原点を通り、また母回帰直線との関係から
\hat{β}_0-β_0=\bar{ε}-\bar{x}(\hat{β_1}-β_1)
従って、分散は次の通り
Var(\hat{β}_0-β_0)=Var(\hat{β}_0)=Var(\bar{ε}-\bar{x}(\hat{β_1}-β_1))
分散の基本公式から
Var(\hat{β}_0)=Var(\bar{ε})+\bar{x}^2Var(\hat{β}_1)−2\bar{x}Cov(\bar{ε},\hat{β}_1−β_1)
共分散は計算すると0となるため、
Var(\hat{β}_0)=\frac{σ^2}{n}+\bar{x}^2*\frac{σ^2}{T_{xx}}=\boxed{σ^2(\frac{1}{n}+\frac{\bar{x}^2}{T_{xx}})}
- 回帰直線におけるyの信頼区間
y_i=β_0+β_1x_i+ε_i
予測値は以下の通りなので、分散はこれまでの議論と共分散を用いて次のように計算できる。
\hat{y_i}=\hat{β}_0+\hat{β}_1x_i
V(\hat{y_i})=Var(\hat{β}_0)+x_i^2Var(\hat{β}_1)+2x_iCov(\hat{β}_0,\hat{β}_1)=σ^2(\frac{1}{n}+\frac{\bar{x}^2}{T_{xx}})+\frac{σ^2}{T_{xx}}-2x_i\bar{x}\frac{σ^2}{T_{xx}}
=σ^2\Bigl(\frac{1}{n}+\frac{(x_i-\bar{x})^2}{T_{xx}}\Bigr)
以上から$y$の標準誤差が与えられるため、母数$y$の信頼区間を次のように求められる。ここで未知$σ$を$σ \rightarrow \hat{σ}$と置き換える。
\boxed{y=\hat{y}±t_{n-2,0.025}×\hat{σ}\sqrt{\frac{1}{n}+\frac{(x_i-\bar{x})^2}{T_{xx}}}}
不偏分散$\hat{σ}^2$は自由度$n-2$の推定量であり、以下の通り(再掲)。
\hat{σ}^2=\frac{\sum_{i=1}^n(\hat{y_i}-y_i)^2}{n-2}=\frac{T_{yy}}{n-2}
5.1.3 回帰係数に関する検定
単回帰モデルでは、係数$β$について帰無仮説$H_0:β=β_0(=0)$が用いられる。説明変数$x$が目的変数$y$の変化に影響を与えない、という帰無仮説と言い換えられる。検証方法は単純で$5.1.2のt$統計検定量を求めて、$t_{0.025}(n-2)$と比較すればよいだけであり、有意水準5%であれば以下の式で帰無仮説を棄却できる。
t=\frac{\hat{β}}{σ/\sqrt{T_{xx}} } ≧ t_{0.025}(n-2)
5.1.4 平均への回帰
繰り返し観測すると、2回目は1回目より全体の平均へ近づく傾向があることを、平均への回帰と呼ぶ。回帰分析が持っているこの性質を見落として、1回目と2回目のアクションの差によるものと誤って判断することを、回帰の錯誤と呼ぶ。
5.1.5 線形重回帰モデル
線形重回帰モデルは次の通り。$β_i$は偏回帰係数と呼ぶ。偏回帰係数がその説明変数を1単位増加させたときに、それ以外の説明変数の値を固定した場合に、$y$がどの程度増加するかを示していると解釈できる。しかし、実際には説明変数間に相関があるので、ある説明変数だけを変化させようとしても、他の説明変数も変化する。
y_i=β_0+β_1x_{1i}+\cdots +β_px_{pi}+ε_i
回帰係数の推定方法は単回帰同様に最小二乗法であり、目的変数の実測値と予測値の差を最小にする仮定をする。
S(\hat{\beta}_0, \hat{\beta}_1, \cdots ,\hat{\beta}_p ) = \sum_{i=1}^n \left( y_i - (\hat{\beta}_0 + \cdots + \hat{\beta}_px_{pi}) \right)^2
説明変数$\hat{β}_i$ごとに偏微分して0とおいた式(正規方程式)を求めて整理すると、$p$元連立1次方程式が得られる。
\hat{σ}_{11}\hat{β}_1+ \cdots +\hat{σ}_{1p}\hat{β}_p=\hat{σ}_{1y}
\vdots
\hat{σ}_{p1}\hat{β}_1+ \cdots +\hat{σ}_{pp}\hat{β}_p=\hat{σ}_{py}
および
\hat{β}_0=\bar{y}-(\hat{\beta}_1\bar{x}_1 + \cdots + \hat{\beta}_p\bar{x}_{pi})
これらから、偏回帰係数を求めて、最後に切片を計算する。但し、$\hat{σ}_{ij}$は$x_i$と$x_j$の共分散の
不偏推定量、$\hat{σ}_{iy}$は$x_i$と$y$の共分散の不偏推定量である。また、期待値は以下の通りである。各偏回帰係数は不偏推定量である。
E[\hat{β}_j]=β_j, \quad j=0,1,\cdots, p
E[\hat{y}_i]=E[\hat{β}_0+ \cdots +\hat{β}_px_{pi}]=β_0+β_1x_{1i}+\cdots +β_px_{pi}
また単回帰同様に重回帰にも以下の性質が成り立つ。
- 予測値の平均と観測値の平均は等しい
- 残差$e_i=y_i-\hat{y}_i$の平均は0である
- 下の回帰式は点$(\bar{x}_1,+ \cdots +\bar{x}_p,\bar{y})$を通る
y=\hat{β}_0+\hat{β}_1x_{1}+ \cdots +\hat{β}_px_{p}
- 予測値$\hat{y}$と残差$e_i=y_i-\hat{y}_i$の相関係数は$0$である
- 平方和は次の通り分解できる
\sum_i(y_i-\bar{y})^2=\sum_i(\hat{y}_i-\bar{y})^2+\sum_i(y_i-\hat{y}_i)^2
行列形式の重回帰の表現
1. 行列とベクトルの基礎
- 次元（形）：$A \in \mathbb{R}^{m \times n}$ は $m$ 行 $n$ 列。積 $AB$ が定義されるには、左の列数 = 右の行数。
- 転置の性質：$(AB)^{\top} = B^{\top} A^{\top}$。
- 逆行列：正方行列 $A$ が可逆なら $A^{-1}$ が存在し、$AA^{-1} = I$。
- 階数（rank）：列ベクトルの独立本数。
- （半）正定値：$A$ が対称で、任意の $v \neq 0$ で $v^\top A v > 0$（$\ge 0$ なら半正定値）。$X^\top X$ は常に半正定値。
- ノルム：$\lVert v \rVert^2 = v^\top v$
2. 重回帰の行列表現
データ数 $n$、説明変数数 $p$。切片を最左列に入れた設計行列を
X =
\begin{bmatrix}
1 & x_{11} & \dots & x_{1p} \\
\vdots & \vdots & & \vdots \\
1 & x_{n1} & \dots & x_{np}
\end{bmatrix}
\in \mathbb{R}^{n \times (p+1)},
\quad
y =
\begin{bmatrix}
y_1 \\ \vdots \\ y_n
\end{bmatrix}
モデルは
\boxed{y = X\beta + \varepsilon}
- $\beta = [\beta_0,\beta_1,\dots,\beta_p]^\top$（$\beta_0$ は切片）
- 仮定：$\mathbb{E}[\varepsilon]=0,\ \mathrm{Var}(\varepsilon)=\sigma^2 I$
3. 最小二乗法
目的：残差平方和 $S(\beta)$ を最小化
S(\beta) = \lVert y - X\beta \rVert^2
= (y - X\beta)^\top (y - X\beta)
- 3.1 展開
スカラーの性質$(a^\top b = b^\top a$）と $(AB)^\top = B^\top A^\top$ を用いて
\begin{aligned}
S(\beta)
= y^\top y - y^\top X\beta - \beta^\top X^\top y + \beta^\top X^\top X \beta \\
= y^\top y - 2\,\beta^\top X^\top y + \beta^\top X^\top X \beta
\end{aligned}
- 3.2 一次条件（勾配 = 0）
行列微分の基本：$\frac{\partial}{\partial \beta}(\beta^\top A \beta) = (A + A^\top)\beta$、
$\frac{\partial}{\partial \beta}(c^\top \beta) = c$。ここで $X^\top X$ は対称だから
\frac{\partial S}{\partial \beta}
= -2 X^\top y + 2 X^\top X \beta
最小点では $\frac{\partial S}{\partial \beta} = 0$ なので
\boxed{ X^\top X \hat{\beta} = X^\top y } \quad \text{（正規方程式）}
- 3.3 二次条件（凸性と一意性）
2階偏微分を与えるヘッセ行列は
\nabla_\beta^2 S(\beta) = 2 X^\top X
$X^\top X$ は半正定値（任意の $v$ に対し $v^\top X^\top X v = \lVert Xv\rVert^2 \ge 0$）。
$X$ が フルランク なら $X^\top X$ は 正定値 で可逆、よって以下の通り一意解を得る。
\boxed{\hat{\beta} = (X^\top X)^{-1} X^\top y}
（ランク不足の場合、最小解は一意でない。ノルム最小解は擬似逆 $X^+$ を使って $\hat{\beta} = X^+ y$）
4. 推定量の分布と分散
モデル $y = X\beta + \varepsilon$ に対し、
\hat{\beta}
= (X^\top X)^{-1} X^\top y
= \beta + (X^\top X)^{-1} X^\top \varepsilon
もし $\mathbb{E}[\varepsilon] = 0,\ \mathrm{Var}(\varepsilon) = \sigma^2 I$ なら
\boxed{ \mathbb{E}[\hat{\beta}] = \beta, \quad
\mathrm{Var}(\hat{\beta}) = \sigma^2 (X^\top X)^{-1} }
（不偏かつ分散は $(X^\top X)^{-1}$ に比例。条件数が悪いと係数が不安定。）
- 残差平方和 $RSS = \lVert y - \hat{y}\rVert^2$
- 不偏推定量 $\hat{\sigma}^2 = RSS/(n - p - 1)$（ここで $p$ は説明変数の数、切片を含めるとパラメータは $p+1$）
- 新規点 $x$ の点予測 $\widehat{y} =$ $x^\top \hat{\beta}$ の分散は
\mathrm{Var}(\widehat{y}) = \sigma^2 x^\top (X^\top X)^{-1} x
- 分散の導出
推定誤差を明示する
\hat{\beta} - \beta
= (X^\top X)^{-1} X^\top y - \beta
= (X^\top X)^{-1} X^\top (X\beta + \varepsilon) - \beta
= (X^\top X)^{-1} X^\top \varepsilon
したがって、
\begin{aligned}
\operatorname{Var}(\hat{\beta})
= \operatorname{Var}\bigl((X^\top X)^{-1} X^\top \varepsilon\bigr)
\end{aligned}
行列の分散の性質：$\operatorname{Var}(A Z) = A \operatorname{Var}(Z) A^\top$（定数行列 $A$ と確率ベクトル $Z$）を利用する(証明は後述)。ここで $A = (X^\top X)^{-1} X^\top$、$Z = \varepsilon$であるから
\operatorname{Var}(\hat{\beta})
= (X^\top X)^{-1} X^\top \operatorname{Var}(\varepsilon) X (X^\top X)^{-1}
= (X^\top X)^{-1} X^\top (\sigma^2 I) X (X^\top X)^{-1}
= \sigma^2 (X^\top X)^{-1} X^\top X (X^\top X)^{-1}
= \boxed{\sigma^2 (X^\top X)^{-1}}
- "行列の分散の性質"の導出
確率ベクトル $Z$ と定数行列 $A$ に対して、線形変換 $AZ$ の分散が
\operatorname{Var}(AZ) = A\operatorname{Var}(Z)A^\top
となることを導出する。
前提
- $Z$：$n \times 1$ の確率ベクトル
- $A$：$m \times n$ の定数行列
- $AZ$：$m \times 1$ の新しい確率ベクトル
ベクトルの分散共分散行列は次で定義される
\operatorname{Var}(Z)
= \mathbb{E}\Big[(Z - \mathbb{E}[Z])(Z - \mathbb{E}[Z])^\top\Big]
$\operatorname{Var}(AZ)$ の定義
- 線形変換 $AZ$ の分散の定義は
\operatorname{Var}(AZ)
= \mathbb{E}\Big[(AZ - \mathbb{E}[AZ])(AZ - \mathbb{E}[AZ])^\top\Big]
期待値 $\mathbb{E}[AZ]$ を整理
期待値の線形性より$\mathbb{E}[AZ] = A\mathbb{E}[Z]$であるので
AZ - \mathbb{E}[AZ] = A(Z - \mathbb{E}[Z])
分散の中身を展開する
定義に代入する
\operatorname{Var}(AZ)
= \mathbb{E}\Big[ A(Z - \mathbb{E}[Z]) \cdot (A(Z - \mathbb{E}[Z]))^\top \Big]
転置の性質 $(AB)^\top = B^\top A^\top$ を使うと
(A(Z - \mathbb{E}[Z]))^\top
= (Z - \mathbb{E}[Z])^\top A^\top
よって
\operatorname{Var}(AZ)
= \mathbb{E}\Big[A (Z - \mathbb{E}[Z])(Z - \mathbb{E}[Z])^\top A^\top\Big]
定数行列は期待値の外に出す
$A$ と $A^\top$ は定数行列であり確率変数ではないため、期待値の外に出せる
\operatorname{Var}(AZ)
= A\mathbb{E}\Big[(Z - \mathbb{E}[Z])(Z - \mathbb{E}[Z])^\top\Big]A^\top
よって
\boxed{\operatorname{Var}(AZ)
= A\operatorname{Var}(Z)A^\top}
5. 切片・中心化・標準化
- 切片を入れるなら $X$ の先頭列をすべて $1$ にする。
- 各列を 中心化（平均0）すると $\hat{\beta}_0 = \bar{y}$ となり解釈が容易。
- 標準化（分散1）はスケール差の調整・正則化の前処理に有効。
- 多重共線性が強いと $X^\top X$ が悪条件化して分散が増大 → Ridge/Lasso や設計の見直し。
6. 重回帰：平均応答の信頼区間と予測区間
6-1. 新しい説明ベクトル $(x_{0})$ における平均応答の信頼区間
点推定（平均応答）
\hat{\mu}_{0} = x_{0}^{\top}\hat{\beta}
であるから、行列の分散の性質：$\operatorname{Var}(A Z) = A \operatorname{Var}(Z) A^\top$より、標準誤差は
\operatorname{SE}(\widehat{\mu}_{0}) = \hat\sigma\sqrt{x_{0}^{\top}(X^{\top}X)^{-1}x_{0}}
信頼水準$α$、自由度 $(n-p) $の $t$ 分布に基づく信頼区間
\boxed{\widehat{\mu}_{0}\ \pm\ t_{n-p, 1-\alpha/2}\ \hat\sigma\sqrt{x_{0}^{\top}(X^{\top}X)^{-1}x_{0}}}
6-2. 予測区間 観測値 $(y_{0})$
予測区間は観測ノイズが加わるため、平均応答の区間より広くなります。
\boxed{\widehat{\mu}_{0}\ \pm\ t_{n-p, 1-\alpha/2}\ \hat\sigma\sqrt{1 + x_{0}^{\top}(X^{\top}X)^{-1}x_{0}}}
偏回帰係数の有意性の検定
帰無仮説$H_0:偏回帰係数β_i=0$として検定をする。単回帰同様に、nはサンプルサイズを、kは説明変数の数として、偏回帰係数を標準誤差で割った値について、自由度$n-k-1のt$分布を用いて検定をする。すなわち、次の式の通りである。デザイン行列を$X$とすると、$T_{xx}$が$(X'X)^{-1}$に置き換わる。
\frac{\hat{β_i}-0}{\sqrt{\sigma^2\left[(X'X)^{-1}\right]_{jj}}}