重回帰モデル ¶ N N N 個のサンプルがあり、i i i 番目のサンプルについての回帰式を次のように表記する
Y i = X i T β + u i , i = 1 , … , N Y_i = X_i^T \beta + u_i, \hspace{1em} i=1,\dots, N Y i = X i T β + u i , i = 1 , … , N ここでX i ∈ R D X_i \in \mathbb{R}^D X i ∈ R D は説明変数の行列X ∈ R N × D X\in \mathbb{R}^{N\times D} X ∈ R N × D を1行取り出したもので、Y i , u i ∈ R Y_i, u_i \in \mathbb{R} Y i , u i ∈ R も1つのサンプルの被説明変数と誤差項である。
行列表記にすると
Y = X β + u Y = X \beta + u Y = Xβ + u と表すことができる。
OLS推定量 ¶ 目的関数は残差の二乗和であるため、
L ( β ) = u T u = ( Y − X β ) T ( Y − X β ) = Y T Y − Y T X β − ( X β ) T Y + ( X β ) T X β = Y T Y − 2 β T X T Y + β T X T X β \begin{align}
L(\beta) &= u^T u\\
&= (Y - X \beta)^T (Y - X \beta)\\
&= Y^T Y - Y^T X \beta - (X\beta)^T Y + (X \beta)^T X \beta\\
&= Y^T Y - 2 \beta^T X^T Y + \beta^T X^T X \beta\\
\end{align} L ( β ) = u T u = ( Y − Xβ ) T ( Y − Xβ ) = Y T Y − Y T Xβ − ( Xβ ) T Y + ( Xβ ) T Xβ = Y T Y − 2 β T X T Y + β T X T Xβ である。これを微分してゼロとおくと
∂ L ( β ) ∂ β = − 2 X T Y + 2 X T X β = 0 \frac{\partial L(\beta)}
{\partial \beta}
= - 2 X^T Y + 2 X^T X \beta
= 0 ∂ β ∂ L ( β ) = − 2 X T Y + 2 X T Xβ = 0 となり、β \beta β について解くと
β = ( X T X ) − 1 X T Y \beta
= (X^T X)^{-1} X^T Y β = ( X T X ) − 1 X T Y 式展開メモ
ベクトルx x x 、行列A A A に対し
∂ A x ∂ x T = A \frac{\partial A x}{\partial x^T} = A ∂ x T ∂ A x = A ∂ x T B x ∂ x = ( B + B T ) x \frac{\partial x^T B x}{\partial x} = (B + B^T) x ∂ x ∂ x T B x = ( B + B T ) x もしBが対称行列なら( B + B T ) x = 2 B x (B + B^T) x = 2Bx ( B + B T ) x = 2 B x
import numpy as np
x0 = np.array([1, 1, 1])
x1 = np.array([1, 2, 3])
x2 = np.array([2, 8, 9])
X = np.array([x0, x1, x2]).T
beta = np.array([3, 5, 7]) # 真のbeta
y = X @ beta
# OLS推定量
beta_hat = np.linalg.inv(X.T @ X) @ (X.T @ y)
beta_hat.round(1)OLS推定量の別表記 ¶ 一致性や不偏性の議論のための準備として、OLS推定量を変形する。
重回帰モデルY = X β + u Y = X\beta + u Y = Xβ + u をOLS推定量β ^ = ( X T X ) − 1 X T Y \hat{\beta} = (X^T X)^{-1} X^T Y β ^ = ( X T X ) − 1 X T Y に代入して変形すると
β ^ = ( X T X ) − 1 X T Y = ( X T X ) − 1 X T ( X β + u ) = ( X T X ) − 1 X T X β + ( X T X ) − 1 X T u = β + ( X T X ) − 1 X T u = β + ( 1 N X T X ) − 1 1 N X T u \begin{align}
\hat{\beta}
&= (X^T X)^{-1} X^T Y\\
&= (X^T X)^{-1} X^T (X\beta + u)\\
&= (X^T X)^{-1} X^T X\beta + (X^T X)^{-1} X^T u\\
&= \beta + (X^T X)^{-1} X^T u\\
&= \beta + \left(\frac{1}{N} X^T X\right)^{-1} \frac{1}{N} X^T u
\end{align} β ^ = ( X T X ) − 1 X T Y = ( X T X ) − 1 X T ( Xβ + u ) = ( X T X ) − 1 X T Xβ + ( X T X ) − 1 X T u = β + ( X T X ) − 1 X T u = β + ( N 1 X T X ) − 1 N 1 X T u となる。
単回帰でいうと
β ^ = β + C o v ( X , u ) V a r ( X ) \hat{\beta} = \beta + \frac{Cov(X, u)}{Var(X)} β ^ = β + Va r ( X ) C o v ( X , u ) である。
単回帰モデル Y = α + β X + u Y=\alpha+\beta X+u Y = α + βX + u は両辺をX X X と共分散をとると
Cov ( X , Y ) = Cov ( X , α + β X + u ) = Cov ( X , α ) ⏟ 定数との共分散 → 0 + β Cov ( X , X ) ⏟ = Var ( X ) + Cov ( X , u ) = β Var ( X ) + Cov ( X , u ) \begin{aligned}
\operatorname{Cov}(X, Y)
&= \operatorname{Cov}(X, \alpha + \beta X + u)\\
&= \underbrace{ \operatorname{Cov}(X, \alpha) }_{定数との共分散 → 0}
+ \beta \underbrace{ \operatorname{Cov}(X, X) }_{=\operatorname{Var}(X) }
+\operatorname{Cov}(X, u) \\
&= \beta \operatorname{Var}(X) + \operatorname{Cov}(X, u)\\
\end{aligned} Cov ( X , Y ) = Cov ( X , α + βX + u ) = 定数との共分散 → 0 Cov ( X , α ) + β = Var ( X ) Cov ( X , X ) + Cov ( X , u ) = β Var ( X ) + Cov ( X , u ) 傾き係数のOLS推定量に代入すると
β ^ = Cov ( X , Y ) Var ( X ) = β Var ( X ) + Cov ( X , u ) Var ( X ) = β + Cov ( X , u ) Var ( X ) \begin{aligned}
\hat{\beta} &= \frac{\operatorname{Cov}(X, Y)}{\operatorname{Var}(X)}\\
&= \frac{ \beta \operatorname{Var}(X) + \operatorname{Cov}(X, u) }{\operatorname{Var}(X)}\\
&= \beta + \frac{\operatorname{Cov}(X, u)}{\operatorname{Var}(X)}
\end{aligned} β ^ = Var ( X ) Cov ( X , Y ) = Var ( X ) β Var ( X ) + Cov ( X , u ) = β + Var ( X ) Cov ( X , u ) X X X に内生性がある、すなわちC o v ( X , u ) ≠ 0 Cov(X, u) \neq 0 C o v ( X , u ) = 0 であるとβ ^ ≠ β \hat{\beta} \neq \beta β ^ = β となる
OLS推定量のバリアンス ¶ β ^ = β + ( X T X ) − 1 X T u \hat{\beta} = \beta + (X^T X)^{-1} X^T u β ^ = β + ( X T X ) − 1 X T u は u ∼ N ( 0 , σ 2 I ) u\sim N(0, \sigma^2 I) u ∼ N ( 0 , σ 2 I ) の仮定より、β ^ ∼ N ( β , σ 2 ( X ⊤ X ) − 1 ) \hat{\beta} \sim N(\beta, \sigma^2 (X^\top X)^{-1}) β ^ ∼ N ( β , σ 2 ( X ⊤ X ) − 1 ) となる。
よって Var [ β ^ ] = σ 2 ( X ⊤ X ) − 1 \operatorname{Var}[\hat{\beta}] = \sigma^2 (X^\top X)^{-1} Var [ β ^ ] = σ 2 ( X ⊤ X ) − 1 となる
以下の定理を使う
定理
u u u を確率変数ベクトルとし、μ ∈ R n \mu\in\mathbb{R}^n μ ∈ R n 、b ∈ R p b\in\mathbb{R}^p b ∈ R p 、C ∈ R p × n C\in\mathbb{R}^{p\times n} C ∈ R p × n 、rank C = p \operatorname{rank} C = p rank C = p とする。
u ∼ N ( μ , Σ ) u \sim N(\mu, \Sigma) u ∼ N ( μ , Σ ) のとき、C u + b ∼ N ( C μ + b , C Σ C ⊤ ) Cu + b \sim N(C\mu + b, C\Sigma C^\top) C u + b ∼ N ( C μ + b , C Σ C ⊤ )
C : = ( X T X ) − 1 X T C := (X^T X)^{-1} X^T C := ( X T X ) − 1 X T とおけば
C T = [ ( X T X ) − 1 X T ] T = X [ ( X T X ) − 1 ] T ( ∵ ( A B ) T = B T A T ) = X [ ( X T X ) T ] − 1 ( ∵ ( A − 1 ) T = ( A T ) − 1 ) = X ( X T X ) − 1 ( ∵ ( X T X ) T = X T X ) \begin{aligned}
C^T = [(X^T X)^{-1} X^T]^T
&= X [(X^T X)^{-1}]^T \quad (\because (AB)^T = B^T A^T) \\
&= X [(X^T X)^T]^{-1} \quad (\because (A^{-1})^T = (A^T)^{-1} ) \\
&= X (X^T X)^{-1} \quad (\because (X^T X)^T=X^T X ) \\
\end{aligned} C T = [( X T X ) − 1 X T ] T = X [( X T X ) − 1 ] T ( ∵ ( A B ) T = B T A T ) = X [( X T X ) T ] − 1 ( ∵ ( A − 1 ) T = ( A T ) − 1 ) = X ( X T X ) − 1 ( ∵ ( X T X ) T = X T X ) であるため、u ∼ N ( 0 , σ 2 I ) u\sim N(0, \sigma^2 I) u ∼ N ( 0 , σ 2 I ) の仮定が満たされるとき、
β + C u ∼ N ( β , σ 2 C C T ) = β + ( X T X ) − 1 X T u ∼ N ( β , σ 2 ( X T X ) − 1 X T X ( X T X ) − 1 ) = β + ( X T X ) − 1 X T u ∼ N ( β , σ 2 ( X T X ) − 1 ) \begin{aligned}
&\beta + C u \sim N(\beta, \sigma^2 CC^T)\\
&= \beta + (X^T X)^{-1} X^T u \sim
N(\beta, \sigma^2 (X^T X)^{-1} X^T X (X^T X)^{-1} )\\
&= \beta + (X^T X)^{-1} X^T u \sim
N(\beta, \sigma^2 (X^T X)^{-1})
\end{aligned} β + C u ∼ N ( β , σ 2 C C T ) = β + ( X T X ) − 1 X T u ∼ N ( β , σ 2 ( X T X ) − 1 X T X ( X T X ) − 1 ) = β + ( X T X ) − 1 X T u ∼ N ( β , σ 2 ( X T X ) − 1 ) i.i.d.:( Y , X ) (Y, X) ( Y , X ) は独立同一分布に従う
外生性:E [ u ∣ X ] = 0 E[u|X] = 0 E [ u ∣ X ] = 0
異常値がない:X , u X, u X , u は4次までのモーメントを持つ
多重共線性がない:任意の∑ j = 0 k a j 2 = 1 \sum_{j=0}^k a_j^2 = 1 ∑ j = 0 k a j 2 = 1 となるa 0 , … , a k a_0,\dots,a_k a 0 , … , a k についてE [ ( a 0 + a 1 X 1 + ⋯ + a k X k ) 2 ] > 0 E[(a_0 + a_1 X_1 + \cdots + a_k X_k)^2]>0 E [( a 0 + a 1 X 1 + ⋯ + a k X k ) 2 ] > 0 が成り立つ
不偏性 ¶ 外生性 ¶ 単回帰モデルY = α + β X + u Y = \alpha + \beta X + u Y = α + βX + u を例にとる。
説明変数X X X と誤差項u u u が
E ( u ∣ X ) = 0 E(u | X) = 0 E ( u ∣ X ) = 0 を満たすとき、X X X は外生変数であるという。
また、外生性の条件は別の表現もできる
説明変数X X X が外生変数ならば、
E ( u ) = 0 , E ( X , u ) = 0 E(u) = 0, E(X, u) = 0 E ( u ) = 0 , E ( X , u ) = 0 E ( X u ) = E X [ E ( X u ∣ X ) ] = E X [ X E ( u ∣ X ) ⏟ = 0 ] = E X ( X ⋅ 0 ) = 0 E ( u ) = E X [ E ( u ∣ X ) ⏟ = 0 ] = 0 E(X u) = E_{X}[ E(X u | X) ] = E_{X}[ X \underbrace{ E(u|X) }_{ =0 } ] = E_{X}(X\cdot 0) = 0
\\
E(u) = E_{X}[ \underbrace{ E(u|X) }_{ =0 } ] = 0 E ( X u ) = E X [ E ( X u ∣ X )] = E X [ X = 0 E ( u ∣ X ) ] = E X ( X ⋅ 0 ) = 0 E ( u ) = E X [ = 0 E ( u ∣ X ) ] = 0 さらに、共分散との関係も導出できる
説明変数X X X が外生変数ならば、
C o v ( X , u ) = 0 Cov(X, u) = 0 C o v ( X , u ) = 0 C o v ( X , u ) = E ( X u ) ⏟ = 0 − E ( X ) E ( u ) ⏟ = 0 = 0 Cov(X, u) = \underbrace{ E(X u) }_{ =0 }
- E(X) \underbrace{ E(u) }_{ =0 }
= 0 C o v ( X , u ) = = 0 E ( X u ) − E ( X ) = 0 E ( u ) = 0 OLS推定量の不偏性 ¶ 単回帰モデルY = α + β X + u Y = \alpha + \beta X + u Y = α + βX + u のOLS推定量β \beta β の確率極限は
Undefined control sequence: \DeclareMathOperator at position 1: \̲D̲e̲c̲l̲a̲r̲e̲M̲a̲t̲h̲O̲p̲e̲r̲a̲t̲o̲r̲{\cov}{\text{Co…
\DeclareMathOperator{\cov}{\text{Cov}}
\DeclareMathOperator{\var}{\text{Var}}
\DeclareMathOperator{\plim}{\text{plim}}
\plim \beta = \beta + \frac{\cov(X, u)}{\var(X)}となる。外生性が満たされるとき \cov(X, u) = 0 であるため、plim β = β \plim \beta = \beta plim β = β となり、OLS推定量は母回帰係数の不偏推定量となる。
Undefined control sequence: \DeclareMathOperator at position 1: \̲D̲e̲c̲l̲a̲r̲e̲M̲a̲t̲h̲O̲p̲e̲r̲a̲t̲o̲r̲{\cov}{\text{Co…
\DeclareMathOperator{\cov}{\text{Cov}}
\DeclareMathOperator{\var}{\text{Var}}
\DeclareMathOperator{\plim}{\text{plim}}
\begin{align}
\plim \beta
&= \frac{\cov(X, Y)}{\var(X)}\\
&= \frac{\cov(X, \alpha + \beta X + u)}{\var(X)}\\
&= \frac{\cov(X, \alpha) + \cov(X, \beta X) + \cov(X, u) }{\var(X)}\\
&= \frac{\beta \var(X) + \cov(X, u) }{\var(X)}\\
&= \beta + \frac{\cov(X, u)}{\var(X)}
\end{align}標本値、確率変数の和は、加える前の個々の共分散の和になる。すなわち、共分散においては分配法則が成り立つ。
Undefined control sequence: \cov at position 1: \̲c̲o̲v̲(X + Z, Y) = \c…
\cov(X + Z, Y) = \cov(X, Y) + \cov(Z, Y)参考:確率統計 – 分散と共分散 – TauStation
α , β \alpha, \beta α , β は定数と考えると、
Undefined control sequence: \cov at position 1: \̲c̲o̲v̲(X, \alpha)
= E…
\cov(X, \alpha)
= E(X \alpha) - E(X) E(\alpha)\\
= \alpha E(X) - \alpha E(X)\\
= 0であり
Undefined control sequence: \cov at position 1: \̲c̲o̲v̲(X, \beta X)
= …
\cov(X, \beta X)
= E(X \beta X) - E(X) E(\beta X)\\
= \beta E(X^2) - \beta E(X)^2\\
= \beta \var(X)一致性 ¶ 異常値がない(X , u X, u X , u は4次までのモーメントを持つ)という仮定と大数の法則により以下が成立する
1 N X T X = 1 N ∑ i = 1 N X i X i T ⟶ p E ( X i X i T ) 1 N X T u = 1 N ∑ i = 1 N X i u i ⟶ p E ( X i u i ) = 0 \begin{align}
\frac{1}{N} X^T X
&= \frac{1}{N} \sum^N_{i=1} X_i X_i^T
\overset{p}{\longrightarrow}
E(X_i X_i^T)
\\
\frac{1}{N} X^T u
&= \frac{1}{N} \sum^N_{i=1} X_i u_i
\overset{p}{\longrightarrow}
E(X_i u_i)
= 0
\end{align} N 1 X T X N 1 X T u = N 1 i = 1 ∑ N X i X i T ⟶ p E ( X i X i T ) = N 1 i = 1 ∑ N X i u i ⟶ p E ( X i u i ) = 0 多重共線性がないという仮定により( E ( X i X i T ) ) − 1 (E(X_i X_i^T))^{-1} ( E ( X i X i T ) ) − 1 が存在する
( 1 N X T X ) − 1 1 N X T u ⟶ p 0 \left(
\frac{1}{N} X^T X
\right)^{-1}
\frac{1}{N} X^T u
\overset{p}{\longrightarrow}
0 ( N 1 X T X ) − 1 N 1 X T u ⟶ p 0 よって
β ^ ⟶ p β \hat{\beta}
\overset{p}{\longrightarrow}
\beta β ^ ⟶ p β 漸近正規性 ¶ OLS推定量
β ^ = β + ( 1 N X T X ) − 1 1 N X T u \hat{\beta}
= \beta + \left(\frac{1}{N} X^T X\right)^{-1} \frac{1}{N} X^T u β ^ = β + ( N 1 X T X ) − 1 N 1 X T u を整理して以下の形にする
N ( β ^ − β ) = ( 1 N X T X ) − 1 1 N X T u \sqrt{N} (\hat{\beta} - \beta)
=
\left(
\frac{1}{N} X^T X
\right)^{-1}
\frac{1}{\sqrt{N}} X^T u N ( β ^ − β ) = ( N 1 X T X ) − 1 N 1 X T u 1 N X T u \frac{1}{\sqrt{N}} X^T u N 1 X T u は1 N ∑ i = 1 N X i u i \frac{1}{\sqrt{N}} \sum^N_{i=1} X_i u_i N 1 ∑ i = 1 N X i u i と書くことができる。OLSの仮定より
E ( X i u i ) = 0 V a r ( X i u i ) = E ( u 2 X i X i T ) \begin{align}
E(X_i u_i) &= 0\\
Var(X_i u_i) &= E(u^2 X_i X^T_i)
\end{align} E ( X i u i ) Va r ( X i u i ) = 0 = E ( u 2 X i X i T ) なので、中心極限定理により
1 N X T u = 1 N ∑ i = 1 N X i u i ⟶ d N ( 0 , E ( u i 2 X i X i T ) ) \frac{1}{\sqrt{N}} X^T u
= \frac{1}{\sqrt{N}} \sum^N_{i=1} X_i u_i
\overset{d}{\longrightarrow} N\left( 0, E(u_i^2 X_i X_i^T) \right) N 1 X T u = N 1 i = 1 ∑ N X i u i ⟶ d N ( 0 , E ( u i 2 X i X i T ) ) となる。
一致性のときに導出した
1 N X T X ⟶ p E ( X i X i T ) \frac{1}{N} X^T X
\overset{p}{\longrightarrow} E(X_i X_i^T) N 1 X T X ⟶ p E ( X i X i T ) を使うと、スルツキーの定理を用いて
N ( β ^ − β ) = ( 1 N X T X ) − 1 1 N X T u ⟶ d ( E ( X i X i T ) ) − 1 × N ( 0 , E ( u i 2 X i X i T ) ) = N ( 0 , ( E ( X i X i T ) ) − 1 E ( u i 2 X i X i T ) ( E ( X i X i T ) ) − 1 ) = N ( 0 , V ) \begin{align}
\sqrt{N} (\hat{\beta} - \beta)
&=
\left(
\frac{1}{N} X^T X
\right)^{-1}
\frac{1}{\sqrt{N}} X^T u
\overset{d}{\longrightarrow} \left( E(X_i X_i^T) \right)^{-1} \times N\left( 0, E(u_i^2 X_i X_i^T) \right)\\
&= N\left( 0,
\left( E(X_i X_i^T) \right)^{-1} E(u_i^2 X_i X_i^T) \left( E(X_i X_i^T) \right)^{-1}
\right)\\
&= N(0, V)\\
\end{align} N ( β ^ − β ) = ( N 1 X T X ) − 1 N 1 X T u ⟶ d ( E ( X i X i T ) ) − 1 × N ( 0 , E ( u i 2 X i X i T ) ) = N ( 0 , ( E ( X i X i T ) ) − 1 E ( u i 2 X i X i T ) ( E ( X i X i T ) ) − 1 ) = N ( 0 , V ) となる。
確率変数の行列Y N , Y , X N , X ∈ R N × N Y_N, Y, X_N, X \in \mathbb{R}^{N\times N} Y N , Y , X N , X ∈ R N × N 、正則行列C ∈ R N × N C \in \mathbb{R}^{N\times N} C ∈ R N × N があるとする。
N → ∞ N\to \infty N → ∞ のとき
X N → d X Y N → d C X_N \xrightarrow{d} X\\
Y_N \xrightarrow{d} C\\ X N d X Y N d C とする。
このとき、以下の結果が成り立ち、これを スルツキーの定理 という
X N + Y N → d X + C X_N + Y_N \xrightarrow{d} X + C X N + Y N d X + C
Y N X N → d C X Y_N X_N \xrightarrow{d} C X Y N X N d CX
Y N − 1 X N → d C − 1 X Y_N^{-1}X_N \xrightarrow{d} C^{-1} X Y N − 1 X N d C − 1 X
V V V は以下のように一致推定できる
V ^ = [ 1 N ∑ i = 1 N X i X i T ] − 1 1 N ∑ i = 1 N u ^ i 2 X i X i T [ 1 N ∑ i = 1 N X i X i T ] − 1 = ( 1 N X T X ) − 1 1 N X T U ^ X ( 1 N X T X ) − 1 \begin{aligned}
\hat{V} & =\left[\frac{1}{N} \sum_{i=1}^N X_i X_i^T\right]^{-1} \frac{1}{N} \sum_{i=1}^N \hat{u}_i^2 X_i X_i^T\left[\frac{1}{N} \sum_{i=1}^N X_i X_i^T\right]^{-1} \\
& =\left(\frac{1}{N} X^T X\right)^{-1} \frac{1}{N} X^T \hat{U} X\left(\frac{1}{N} X^T X\right)^{-1}
\end{aligned} V ^ = [ N 1 i = 1 ∑ N X i X i T ] − 1 N 1 i = 1 ∑ N u ^ i 2 X i X i T [ N 1 i = 1 ∑ N X i X i T ] − 1 = ( N 1 X T X ) − 1 N 1 X T U ^ X ( N 1 X T X ) − 1 ただし、U ^ \hat{U} U ^ は対角要素にu ^ 1 2 , … , u ^ N 2 \hat{u}_1^2, \dots, \hat{u}_N^2 u ^ 1 2 , … , u ^ N 2 を並べた対角行列である。
バイアスとバリアンス ¶ 最小二乗推定量はすべての線形不偏推定量の中で最もバリアンスが小さい(最良である)ことを示すガウス・マルコフの定理というものがある。
バイアスとバリアンス
統計学・機械学習の分野では誤差をバイアス(bias)とバリアンス(variance)に分けて考えることが多い。
バイアスは真の値θ \theta θ とサンプルを変えて推定を繰り返したときの個々の推定結果の平均E ( θ ^ ) E(\hat{\theta}) E ( θ ^ ) との差で、
バリアンスは推定の分散V ( θ ^ ) V(\hat{\theta}) V ( θ ^ ) である。
B i a s = θ − E ( θ ^ ) V a r i a n c e = V ( θ ^ ) Bias = \theta - E(\hat{\theta})\\
Variance = V(\hat{\theta}) B ia s = θ − E ( θ ^ ) Va r ian ce = V ( θ ^ ) 不偏性 ¶ 任意のパラメータの線形結合θ = α ⊤ β \theta=\boldsymbol{\alpha}^\top {\boldsymbol{\beta}} θ = α ⊤ β を考える。例えばf ( x 0 ) = x 0 ⊤ β f(x_0)=x_0^\top \beta f ( x 0 ) = x 0 ⊤ β がこの形である。
この最小二乗推定値は
θ ^ = α ⊤ β ^ = α ⊤ ( X ⊤ X ) − 1 X ⊤ y \hat{\theta} = \boldsymbol{\alpha}^\top \hat{\boldsymbol{\beta}}
= \boldsymbol{\alpha}^\top (\boldsymbol{X}^\top \boldsymbol{X})^{-1} \boldsymbol{X}^\top \boldsymbol{y} θ ^ = α ⊤ β ^ = α ⊤ ( X ⊤ X ) − 1 X ⊤ y で、期待値をとると
E [ θ ^ ] = E [ α ⊤ β ^ ] = E [ α ⊤ ( X ⊤ X ) − 1 X ⊤ y ] = E [ α ⊤ ( X ⊤ X ) − 1 X ⊤ ( X β + ε ) ] = α ⊤ ( X ⊤ X ) − 1 X ⊤ X β + α ⊤ ( X ⊤ X ) − 1 X ⊤ E [ ε ] = α ⊤ ( X ⊤ X ) − 1 X ⊤ X β = α ⊤ β \begin{align}
E[\hat{\theta}]
&= E[\boldsymbol{\alpha}^\top \hat{\boldsymbol{\beta}}]\\
&= E[\boldsymbol{\alpha}^\top (\boldsymbol{X}^\top \boldsymbol{X})^{-1} \boldsymbol{X}^\top \boldsymbol{y}]\\
% &= \boldsymbol{\alpha}^\top (\boldsymbol{X}^\top \boldsymbol{X})^{-1} \boldsymbol{X}^\top E[\boldsymbol{y}]\\
&= E[\boldsymbol{\alpha^\top (X^\top X)^{-1} X^\top (X\beta + \varepsilon) }] \\
&= \boldsymbol{\alpha^\top (X^\top X)^{-1} X^\top X \beta + \alpha^\top (X^\top X)^{-1} X^\top} E[ \varepsilon ] \\
&= \boldsymbol{\alpha^\top (X^\top X)^{-1} X^\top X \beta}\\
&= \boldsymbol{\alpha}^\top \boldsymbol{\beta}
\end{align} E [ θ ^ ] = E [ α ⊤ β ^ ] = E [ α ⊤ ( X ⊤ X ) − 1 X ⊤ y ] = E [ α ⊤ ( X ⊤ X ) − 1 X ⊤ ( Xβ + ε ) ] = α ⊤ ( X ⊤ X ) − 1 X ⊤ Xβ + α ⊤ ( X ⊤ X ) − 1 X ⊤ E [ ε ] = α ⊤ ( X ⊤ X ) − 1 X ⊤ Xβ = α ⊤ β ※誤差項は仮定によりE [ ε ] = 0 E[\boldsymbol{\varepsilon}]=0 E [ ε ] = 0 であるため
となり(※)、θ ^ \hat{\theta} θ ^ が不偏推定量である(E [ θ ^ ] = θ E[\hat{\theta}] = \theta E [ θ ^ ] = θ )ことがわかる。
バリアンス ¶ α ⊤ β \boldsymbol{\alpha}^\top \boldsymbol{\beta} α ⊤ β に対して不偏のまた別の線形推定量c ⊤ y \boldsymbol{c}^\top \boldsymbol{y} c ⊤ y があるとする。
両者の差を
α ⊤ β ^ − c ⊤ y = [ α ⊤ ( X ⊤ X ) − 1 X − c ⊤ ] y = : d ⊤ y \begin{align}
\boldsymbol{\alpha}^\top \boldsymbol{\hat{\beta}} - \boldsymbol{c}^\top \boldsymbol{y}
&= [ \boldsymbol{\alpha}^\top (\boldsymbol{X}^\top \boldsymbol{X})^{-1} \boldsymbol{X} - \boldsymbol{c}^\top ] \boldsymbol{y}\\
&=: \boldsymbol{d}^\top \boldsymbol{y}
\end{align} α ⊤ β ^ − c ⊤ y = [ α ⊤ ( X ⊤ X ) − 1 X − c ⊤ ] y =: d ⊤ y とおく。このとき、不偏性E [ c ⊤ y ] = α ⊤ β E[\boldsymbol{c}^\top \boldsymbol{y}] = \boldsymbol{\alpha}^\top \boldsymbol{\beta} E [ c ⊤ y ] = α ⊤ β から両者の差の期待値はゼロになるべきであり、
E [ d ⊤ y ] = d ⊤ X β = 0 E[\boldsymbol{d}^\top \boldsymbol{y}]
= \boldsymbol{d}^\top \boldsymbol{X} \boldsymbol{\beta} = 0 E [ d ⊤ y ] = d ⊤ X β = 0 が任意のβ \boldsymbol{\beta} β について成り立たなければならないため、
d ⊤ X = 0 \boldsymbol{d}^\top \boldsymbol{X} = \boldsymbol{0} d ⊤ X = 0 が成り立つ。
次に、2つの確率変数X , Y X, Y X , Y に対し
V [ X + Y ] = V [ X ] + 2 Cov [ X , Y ] + V [ Y ] V[X+Y] = V[X] + 2 \text{Cov}[X, Y] + V[Y] V [ X + Y ] = V [ X ] + 2 Cov [ X , Y ] + V [ Y ] が成り立つから、c ⊤ y \boldsymbol{c}^\top \boldsymbol{y} c ⊤ y の分散は
V [ c ⊤ y ] = V [ α ⊤ β ^ − d ⊤ y ] = V [ α ⊤ β ^ ] − 2 Cov [ α ⊤ ( X ⊤ X ) − 1 X ⊤ y , d ⊤ y ] + V [ d ⊤ y ] \begin{align}
V[\boldsymbol{c^\top y}] &= V[\boldsymbol{ \alpha^\top \hat{\beta} - d^\top y }]\\
&= V[\boldsymbol{ \alpha^\top \hat{\beta} }]
- 2 \text{Cov} [\boldsymbol{ \alpha^\top (X^\top X)^{-1} X^\top y }, \boldsymbol{ d^\top y} ]
+ V[\boldsymbol{ d^\top y}]
\end{align} V [ c ⊤ y ] = V [ α ⊤ β ^ − d ⊤ y ] = V [ α ⊤ β ^ ] − 2 Cov [ α ⊤ ( X ⊤ X ) − 1 X ⊤ y , d ⊤ y ] + V [ d ⊤ y ] と表すことができる。
Cov ( A , B ) = E [ ( A − E [ A ] ) ( B − E [ B ] ) ⊤ ] \text{Cov}(A, B) = E[(A - E[A])(B - E[B])^\top] Cov ( A , B ) = E [( A − E [ A ]) ( B − E [ B ] ) ⊤ ] Cov ( a ⊤ y , b ⊤ y ) = E [ ( a ⊤ y − E [ a ⊤ y ] ) ( b ⊤ y − E [ b ⊤ y ] ) ⊤ ] \text{Cov}(a^\top y, b^\top y)
= E[(a^\top y - E[a^\top y])(b^\top y - E[b^\top y])^\top] Cov ( a ⊤ y , b ⊤ y ) = E [( a ⊤ y − E [ a ⊤ y ]) ( b ⊤ y − E [ b ⊤ y ] ) ⊤ ] Cov ( a ⊤ y , b ⊤ y ) = E [ ( a ⊤ y − E [ a ⊤ y ] ) ( b ⊤ y − E [ b ⊤ y ] ) ⊤ ] = E [ ( a ⊤ y − α ⊤ β ) ( b ⊤ y − E [ b ⊤ y ] ) ⊤ ] \text{Cov}(a^\top y, b^\top y)
= E[(a^\top y - E[a^\top y])(b^\top y - E[b^\top y])^\top]\\
= E[(a^\top y - \alpha^\top \beta)(b^\top y - E[b^\top y])^\top] Cov ( a ⊤ y , b ⊤ y ) = E [( a ⊤ y − E [ a ⊤ y ]) ( b ⊤ y − E [ b ⊤ y ] ) ⊤ ] = E [( a ⊤ y − α ⊤ β ) ( b ⊤ y − E [ b ⊤ y ] ) ⊤ ] 次に、2つの確率変数
a ⊤ y = ∑ a i y i , b ⊤ y = ∑ b i y i \boldsymbol{a^\top y} = \sum a_i y_i,
\hspace{2em}
\boldsymbol{b^\top y} = \sum b_i y_i a ⊤ y = ∑ a i y i , b ⊤ y = ∑ b i y i の共分散は、誤差項ε \boldsymbol{\varepsilon} ε が無相関・等分散の仮定V [ ε ] = σ 2 I V[\boldsymbol{\varepsilon}] = \sigma^2 \boldsymbol{I} V [ ε ] = σ 2 I を満たすとき、
Cov [ a ⊤ y , b ⊤ y ] = Cov [ a ⊤ ε , b ⊤ ε ] = ∑ a i b i σ 2 = ( a ⊤ b ) σ 2 \text{Cov}[\boldsymbol{a^\top y, b^\top y}]
= \text{Cov}[\boldsymbol{a^\top \varepsilon, b^\top \varepsilon}]
= \sum a_i b_i \sigma^2
= (\boldsymbol{a^\top b}) \sigma^2 Cov [ a ⊤ y , b ⊤ y ] = Cov [ a ⊤ ε , b ⊤ ε ] = ∑ a i b i σ 2 = ( a ⊤ b ) σ 2 となることから
Cov [ α ⊤ ( X ⊤ X ) − 1 X ⊤ y , d ⊤ y ] = α ⊤ ( X ⊤ X ) − 1 X ⊤ d ⋅ σ 2 \text{Cov}[\boldsymbol{ \alpha^\top (X^\top X)^{-1} X^\top y, d^\top y }]
= \boldsymbol{\alpha^\top (X^\top X)^{-1} X^\top d } \cdot \sigma^2 Cov [ α ⊤ ( X ⊤ X ) − 1 X ⊤ y , d ⊤ y ] = α ⊤ ( X ⊤ X ) − 1 X ⊤ d ⋅ σ 2 となり、d ⊤ X = 0 \boldsymbol{d}^\top \boldsymbol{X} = \boldsymbol{0} d ⊤ X = 0 よりこれは0となる。
よって
V [ c ⊤ y ] = V [ α ⊤ β ^ ] + V [ d ⊤ y ] V[\boldsymbol{ c^\top y }] = V[\boldsymbol{ \alpha^\top \hat{\beta} }] + V[\boldsymbol{ d^\top y }] V [ c ⊤ y ] = V [ α ⊤ β ^ ] + V [ d ⊤ y ] が成り立ち、分散は非負なので
V [ c ⊤ y ] ≥ V [ α ⊤ β ^ ] V[\boldsymbol{ c^\top y }] \geq V[\boldsymbol{ \alpha^\top \hat{\beta} }] V [ c ⊤ y ] ≥ V [ α ⊤ β ^ ] を意味する。
よってα ⊤ β ^ \boldsymbol{ \alpha^\top \hat{\beta} } α ⊤ β ^ は最良線形不偏推定量BLUEである。
OLS推定の幾何学的意味 ¶ OLS推定量
β ^ = ( X ⊤ X ) − 1 X ⊤ y \hat{\boldsymbol{\beta}} = (X^\top X)^{-1} X^\top \boldsymbol{y} β ^ = ( X ⊤ X ) − 1 X ⊤ y をy ^ = X β ^ \hat{\boldsymbol{y}} = X \hat{\boldsymbol{\beta}} y ^ = X β ^ に代入すると
y ^ = X ( X ⊤ X ) − 1 X ⊤ ⏟ P y = P y \hat{\boldsymbol{y}} = \underbrace{ X (X^\top X)^{-1} X^\top }_{P} \boldsymbol{y}
= P \boldsymbol{y} y ^ = P X ( X ⊤ X ) − 1 X ⊤ y = P y つまり、ベクトルy \boldsymbol{y} y を行列P = X ( X ⊤ X ) − 1 X ⊤ P = X (X^\top X)^{-1} X^\top P = X ( X ⊤ X ) − 1 X ⊤ で射影したものとみなすことができる。
この行列P P P は対称行列で、P 2 = P P^2=P P 2 = P となる。この2つの性質を満たす行列を射影行列という。
P 2 = P P = ( X ( X ⊤ X ) − 1 X ⊤ ) ( X ( X ⊤ X ) − 1 X ⊤ ) = X ( X ⊤ X ) − 1 ( X ⊤ X ) ( X ⊤ X ) − 1 X ⊤ = X ( X ⊤ X ) − 1 X ⊤ = P \begin{aligned}
P^2 & = PP \\
& =(X(X^{\top} X)^{-1} X^{\top})(X(X^{\top} X)^{-1} X^{\top}) \\
& =X(X^{\top} X)^{-1}(X^{\top} X)(X^{\top} X)^{-1} X^{\top} \\
& =X(X^{\top} X)^{-1} X^{\top}
=P
\end{aligned} P 2 = PP = ( X ( X ⊤ X ) − 1 X ⊤ ) ( X ( X ⊤ X ) − 1 X ⊤ ) = X ( X ⊤ X ) − 1 ( X ⊤ X ) ( X ⊤ X ) − 1 X ⊤ = X ( X ⊤ X ) − 1 X ⊤ = P 射影行列は、X X X の列空間ℑ X \Im X ℑ X にベクトルを正射影するという性質がある。y \boldsymbol{y} y のℑ X \Im X ℑ X への射影がy ^ \hat{\boldsymbol{y}} y ^ で、垂線の足が誤差u \boldsymbol{u} u となる。
よって、最小二乗法はy \boldsymbol{y} y からℑ X \Im X ℑ X への射影を求める操作であると捉えることができる。
OLSとBLUE / BUE ¶ OLSはBLUEかBUEか?
行列A A A によってb = A y b = Ay b = A y のような線形結合で表現される推定量b b b は線形推定量という。OLSもA = ( X ⊤ X ) − 1 X ⊤ A=(X^\top X)^{-1} X^\top A = ( X ⊤ X ) − 1 X ⊤ とおけば同じ形になっていることがわかる。
OLSはガウス・マルコフの定理でBLUE(線形不偏推定量のなかで最良)だと示された。
Hansen (2022) は線形制約は不要で、線形と非線形の両方のすべての不偏推定量の中で最良(BUE)だと主張した。
一方で Pötscher & Preinerstorfer (2022) はHansen (2022)に対する批判を展開した。
Portnoy (2022) は「一般線形モデルの場合、十分に広い分布族内のすべての分布に対して、不偏な推定量は線形でなければならない」と述べている。
Hansen(2022) の「OLSが最良の不偏推定量(BUE)」が正しいとしても、Portnoy (2022)によれば「一般線形モデルの不偏推定量は線形推定量」なので結局BLUE=BUEになり、いずれの主張も間違っていないことになる。
参考:
参考文献 ¶ 東京大学出版会『統計学入門』
東京大学出版会『自然科学の統計学』
Hastie, T., Tibshirani, R., Friedman, J. H., & Friedman, J. H. (2009). The elements of statistical learning: data mining, inference, and prediction (Vol. 2, pp. 1-758). New York: springer.
西山慶彦, 新谷元嗣, 川口大司, & 奥井亮. (2019). 計量経済学. Yūhikaku.
Hansen, B. E. (2022). A Modern Gauss–Markov Theorem. Econometrica , 90 (3), 1283–1294. 10.3982/ecta19255 Pötscher, B. M., & Preinerstorfer, D. (2022). A Modern Gauss-Markov Theorem? Really? 10.48550/ARXIV.2203.01425 Portnoy, S. (2022). Linearity of Unbiased Linear Model Estimators. The American Statistician , 76 (4), 372–375. 10.1080/00031305.2022.2076743