はじめに
回帰分析について、理解した範囲でまとめます。概念としては結果だけをみて「そうだなー」と分かりますが、実際に式を導出しようとするとつまづきました。背後にある知識が脆弱であるのが利用かなと思います。
流れとして、
- パラメータ\betaの最小自乗推定の解
- パラメータ\betaと分散\sigma^2の最尤推定
- パラメータの推定値の分散は最適(by クラメール・ラオの不等式)
- 最適解の予測値と誤差の直交性
- 分散の不偏性
といきたいと思います。たくさん書いているけれど、結果は至ってシンプルだと思います。
回帰分析のパラメータ推定
問題設定
ある観測値yをp個の説明変数で説明しようと思います。誤差は正規分布に従うとします。
y = \beta_0 + \beta_1 x_1 + \cdots + \beta_p x_p + \epsilon \\
\epsilon \sim N(0, \sigma^2)
今、n個の観測があり、それぞれに y_i = (1, \mathbf{x}_i^T)\mathbf{\beta} + \epsilon_i が成り立つとします。各標本化における誤差は独立であり、同一の分散値であるとします。
これを通常はまとめて
\mathbf{y} = \mathbf{X}\bm{\beta} + \mathbf{\epsilon}
と書くようです。
状況として、(\mathbf{x}_i, y_i) が i=1,\cdots,nとn個与えられたとき、それぞれの要素xがどのようにyに関係しているかを調べ、また、未知のxがあったときにyを予測することに使えるようです。その為に、\mathbf{\beta}と\sigma^2について調べてみましょう。
パラメータ\mathbf{\beta}の推定
二乗誤差の最小化
まず誤差を最小化するパラメータを求めてみましょう。誤差を最小化する点においては勾配が０なので、
\nabla_\beta |\mathbf{y} - \mathbf{X}\mathbf{\beta}|^2 = -2 \mathbf{X}^T (\mathbf{y-X\beta}) = \mathbf{0}
が成り立ちます。これからnormal equation を得られます。
\mathbf{X}^T\mathbf{X\beta} = \mathbf{X^Ty}
行列は対称行列です。更に、同じ行がない、すなわち、同じサンプルを含まないことを仮定すると、行列は正定値であることが保証されます。これより、逆行列をかけることで推定値を得ることができます。
\hat{\mathbf{\beta}} = (\mathbf{X}^T\mathbf{X})^{-1} \mathbf{X}^T \mathbf{y}
最小二乗推定量(least squares estimator, LSE)と呼ばれるようです。
最尤推定
誤差が正規分布に従うとします。y\sim N(\mathbf{x}^T\mathbf{\beta},\sigma^2) です。先程求めた最小二乗推定量は最尤推定量でもあります。
まとめて書くと
\mathbf{y} \sim N(\mathbf{X}\mathbf{\beta}, \sigma^2 \mathbf{I}_n)
対数尤度 \ell(\beta, \sigma^2) = \log L(\beta,\sigma^2) は、
\ell(\beta,\sigma^2) = \frac{1}{2\sigma^2} | \mathbf{y} - \mathbf{X}\mathbf{\beta}|^2 - \frac{n}{2}\log(2\pi\sigma^2)
です。これ、普通の正規分布の最尤推定ですよね。何故わざわざ対数尤度の微分の式を書いたかというと、このあとこいつを更に微分してHessianも計算するからです。
\nabla_\beta \log f(\beta,\sigma^2)= - \frac{1}{2\sigma^2} 2X^T (y-X\beta)
が０であるという、先ほど同じ計算で、\hat{\beta} が得られます。また、分散についても例えば \tau=\frac{1}{\sigma^2} と変数変換して微分係数が０になる点を求めることができます。
\hat{\mathbf{\beta}}_\mathrm{ML} = (\mathbf{X^TX})^{-1} X^T \mathbf{y} \\
\hat{\sigma}^2_\mathrm{ML} = \frac{1}{n} (\mathbf{y} - X\hat{\beta}_\mathrm{ML})
(\mathbf{y} - X\hat{\beta}_\mathrm{ML})^T
を得ます。
パラメータ\beta の最小分散の不偏推定量
どちらから話しても良いのですが、まずパラメータ\beta の不偏性について調べてみます。最小自乗推定量の期待値を計算すると、、、
E[\hat{\beta}_\mathrm{ML}] = (\mathbf{X}^T\mathbf{X})^{-1} \mathbf{X}^T E[y] \\
= (\mathbf{X}^T\mathbf{X})^{-1} \mathbf{X}^T (X\beta + E[\epsilon]) = \beta
不偏推定量であることが分かります。以後、\hat{\beta} と書くことにします。
パラメータ推定量のばらつきである共分散行列についても計算してみましょう。Cov(\mathbf{y},\mathbf{y})=\sigma^2 \mathbf{I}_nであるので、
Cov(\hat{\beta},\hat{\beta}) = (\mathbf{X}^T\mathbf{X})^{-1} \mathbf{X}^T Cov(\mathbf{y},\mathbf{y}) \mathbf{X} (\mathbf{X}^T\mathbf{X})^{-1} = \sigma^2 (X^TX)^{-1}
これは実は不偏推定量の中での分散の最小値です。というのも、これはクラメール・ラオの不等式の下限
Var[\hat{\beta}] \geq I(\beta)^{-1}
になっているからです。Fisherの情報行列は、対数尤度の２階微分(second derivative)で計算できます。期待値を計算する部分は出てきませんでした。
I(\beta)_n = -E[\mathbf{\nabla}_\beta^2 \log L(\mathbf{\beta},\sigma^2)] = \frac{1}{\sigma^2} \mathbf{X}^T \mathbf{X}
から確認できました。
予測値と誤差の直交性と決定係数
この最尤推定であり不偏性もある推定値\betaを用いたときの予測値は
\hat{y} = \mathbf{X}\hat{\mathbf{\beta}} = \mathbf{X} (\mathbf{X}^T\mathbf{X})^{-1} \mathbf{X}^T \mathbf{y}
と書けます。ここで、これをベクトルyから\hat{y}への写像と考えることします。
P_X = \mathbf{X} (\mathbf{X}^T\mathbf{X})^{-1} \mathbf{X}^T
これは観測データyの説明変数Xの春空間への射影 と考えることができます。この写像を表す行列をP_Xと置くと、これは射影行列になります。P_X^2=P_Xです。
この時の誤差は
e = Y - X \hat{\beta} = (I-P_X)y
になります。というわけで、二つの直交補空間ができました。
- 説明空間への射影: \mathit{P_X}: \mathbf{y} \rightarrow \hat{\mathbf{y}} は次元は p+1
- 誤差への射影 \mathbf{I}-\mathit{P_X}: \mathbf{y} \rightarrow \mathbf{e} はn-p-1次元です。説明しきれない空間ですね。。。
となり、見通しが良くなります。目的変数がある n 次元の空間を p+1 次元の説明変数の空間と誤差 n-p-1 次元の空間に分けることができていますね。
決定係数
平均を引き算してみます。成分が全て１の n\times n の行列を \mathbf{J}_n と書くことにします。
\mathbf{y} - \frac{1}{n}\mathbf{J}_n y = (\mathbf{I}-\frac{1}{n}\mathbf{J}) (X\hat{\beta}+\epsilon)
整理すると、平均値からずれについて下記となります。
\mathbf{y} - \bar{\mathbf{y}} = (\mathbf{X} - \bar{\mathbf{X}}\mathbf{\beta} + \mathbf{\epsilon}
説明空間と誤差が直交補空間にあることから、観測変数の変動は回帰変動と誤差変動に分解できます。
\sum_{i=1}^n (y_i-\bar{y})^2 = \sum_{i=1}^n ((x_i - \bar{x_i})^T\hat{\beta}_{1:p})^2 + \sum_{i=1}^n e_i^2
これから、モデルが目的変数をどのくらい説明変数で説明できているかを評価するのに、その回帰変動がしめる割合を考えようということになり、
R^2 = \frac{\mathbf{回帰変動}}{\mathbf{総変動}}
を決定係数(coefficient of determination)と言います。"R square"とも呼ぶみたいですね。実際には、変数を増やせば増やすほど説明力が上がるので、自由度を考慮した調整済み決定係数を用いることが多いようです。
分散\sigma^2の不偏推定量
分散の推定量
\hat{\sigma}^2 = \frac{1}{n} |\mathbf{y} - \mathbf{X}\hat{\mathbf{\beta}}|^2
の期待値を計算してみましょう。不偏推定量になっているかな？
誤差については \mathbf{e}=\mathbf{y} - \mathbf{X}\hat{\mathbf{\beta}} = (\mathbf{I} - P_X)y=(I - P_X)\epsilon ですが、期待値の計算に出てくる確率変数は誤差 \epsilon なので、
E[\mathbf{e}^T \mathbf{e}] = E[|(\mathrm{I}-P_X) \mathbf{\epsilon}|^2] = Tr[\sigma^2 (I-PX)^2]
ここで何故いきなりトレース出てくるの？と思われるかもしれませんが（自分はそうでした^^;)、落ち着いて考えると納得できます。ここでは、\epsilon_i は独立であることを仮定しているので i\neq j において E[\epsilon_i \epsilon_j]=0 であり、E[\epsilon_i^2]=\sigma^2 です。なので２次形式\epsilon^T (I-P_X)^2 \epsilon の計算でnon-zeroになるのは \epsilon_i の次元が等しいところだけになるのでトレースになります。
では頑張ってトレースを計算してみましょう。Tr(I_n)=nは良いとして、P_Xがどうなるか調べてみましょう。
Tr(P_X) = Tr(\mathbf{X} (\mathbf{X}^T\mathbf{X})^{-1} \mathbf{X}^T) = Tr(\mathbf{X}^T\mathbf{X} (\mathbf{X}^T\mathbf{X})^{-1}) = Tr(I_{p+1}) = p+1
です。これも一瞬何これ、と思うかもしれませんが、トレースについて Tr(AB) = Tr(BA) が実はA, Bが正方行列でなくても成立することを自分で確認できれば納得できると思います。Projection matrixのwikipediaの説明を読んだら、
For linear models, the trace of the projection matrix is equal to the rank of
\mathbf{X}, which is the number of independent parameters of the linear model.
とまで書かれていました。
というわけで、Tr(I - P_X)=n-(p+1) の示せたので、まとめると、
E[\mathbf{e}^T \mathbf{e}]=(n-p-1)\sigma^2
となりました。このとから E[\hat{\sigma}^2]=\sigma^2 となってくれる推定値は
\hat{\sigma}^2 = \frac{1}{n-p-1}\mathbf{e}^T \mathbf{e}
であることが分かりました。
これらの期待値の計算で確認した不偏推定量から、推定値が従う標本分布が分かりました。
- \hat{\beta} \sim N(\beta,\sigma^2 (X^TX)^{-1})
- \frac{(n-p-1)\hat{\sigma^2}}{\hat{\sigma}^2} \sim \chi^2(n-p-1)
めでたしめでたし。
まとめ
線形モデルのパラメータとして回帰係数と誤差の分散を推定する計算方法を確認しました。
回帰係数の最尤推定値は不偏推定値でもあり、分散は最小であることが情報行列の計算から確認できました。
誤差の分散の不偏推定量を計算することもできました。また、最適パラメータを用いた予測値の計算は、観測値を説明可能な空間への射影と考えることができ、誤差との直交性から決定係数を当てはまりの指標として使えることを紹介しました。
まぁ、先は長そうです。検定、一般化線形モデル、回帰診断などまだお勉強まとめが必要なことがある。が、ここまでは自分では納得できそうなので良しとします。
Discussion