前置き
線形回帰モデルは,以下のように表せます:
$$
\mathbf{y} = \mathbf{X}\boldsymbol{\beta} + \boldsymbol{\varepsilon},
$$
ここで,
$$
\boldsymbol{\varepsilon} \sim \mathrm{N}\left(\mathbf{0}, \sigma^2\mathbf{I}_n\right)
$$
です。すなわち,誤差項は平均 $\mathbf{0}$,分散 $\sigma^2\mathbf{I}_n$ の正規分布に従うとします。
この線形回帰モデルに対しては,
- 尤度を最大化する 最尤推定 (MLE),
- 係数ベクトル $\boldsymbol{\beta}$ に事前分布を課す 最大事後確率推定(MAP:事後分布のモード)
という二つの立場から解くことができます。
興味深いことに,このモデルでは
- MLE による推定量は OLS (Ordinary Least Squares;最小二乗法) と一致し,
- 正規事前を仮定した MAP 推定量は Ridge 回帰の解 と一致します。
以下では,これらの同値性を順に示します。
OLS = MLE
最尤推定では,$\boldsymbol{\beta}$ は未知の固定パラメータとみなし,
事前分布は設定しません。
観測データ $\mathbf{y}$ に対する尤度は
$$
p(\mathbf{y} \mid \boldsymbol{\beta}, \sigma^2)
= (2\pi\sigma^2)^{-n/2}
\exp\left(
-\frac{1}{2\sigma^2}
\lVert\ \mathbf{y}-\mathbf{X}\boldsymbol{\beta}\rVert^2
\right)
$$
と書けます。
定数を除いた対数尤度は
$$
\log p(\mathbf{y} \mid \boldsymbol{\beta}, \sigma^2)
\propto
-\lVert\ \mathbf{y}-\mathbf{X}\boldsymbol{\beta}\rVert^2
$$
となります。OLS推定(誤差を最小化する)のと同じ式になりました。
あとはこれを$\boldsymbol{\beta}$で微分すれば良いです。
よって解は
$$
\hat{\boldsymbol{\beta}}=
(\mathbf{X}^\mathsf{T}\mathbf{X})^{-1}\mathbf{X}^\mathsf{T}\mathbf{y}.
$$
となります。いわゆるOLS estimatorになりました。
Ridge = MAP
ベイズ的なリッジ回帰(MAP 推定)は,次のモデルに基づきます:
$$
\mathbf{y} \mid \boldsymbol{\beta}, \sigma^2
\sim \mathrm{N}\left(\mathbf{X}\boldsymbol{\beta}, , \sigma^2 \mathbf{I}_n\right),
$$
今回は,$\boldsymbol{\beta}$ にも(正規)分布を仮定します:
$$
\boldsymbol{\beta} \mid \tau^2
\sim \mathrm{N}\left(\mathbf{0}, \tau^2 \mathbf{I}_p\right).
$$
MAP推定では,未知パラメータ $\boldsymbol{\beta}$ を確率変数とみなし,
ベイズの定理により事後分布
$p(\boldsymbol{\beta}\mid\mathbf{y},\sigma^2,\tau^2)$ を求めます。
ベイズの定理は
$$
p(\boldsymbol{\beta}\mid\mathbf{y},\sigma^2,\tau^2)
= \frac{
p(\mathbf{y}\mid\boldsymbol{\beta},\sigma^2),
p(\boldsymbol{\beta}\mid\tau^2)
}{
p(\mathbf{y}\mid\sigma^2,\tau^2)
}
\propto
p(\mathbf{y}\mid\boldsymbol{\beta},\sigma^2)
p(\boldsymbol{\beta}\mid\tau^2)
$$
と書けます。ここで,
$p(\mathbf{y}\mid\boldsymbol{\beta},\sigma^2)$ が尤度,
$p(\boldsymbol{\beta}\mid\tau^2)$ が事前分布に対応します。
したがって,指数部分を合わせ対数を取ると(正規化定数を除いて)
$$
\log p(\boldsymbol{\beta}\mid\mathbf{y},\sigma^2,\tau^2)=-\frac{1}{2\sigma^2}
(\mathbf{y}-\mathbf{X}\boldsymbol{\beta})^{\mathsf{T}}
(\mathbf{y}-\mathbf{X}\boldsymbol{\beta})-\frac{1}{2\tau^2}
\boldsymbol{\beta}^{\mathsf{T}}\boldsymbol{\beta}+ \text{const}.
$$
二乗ノルムの表記にすると事後分布の対数は(定数項を除いて)
$$
\log p(\boldsymbol{\beta}\mid\mathbf{y},\sigma^2,\tau^2)
\propto -(\frac{1}{2\sigma^2}
\left\lVert \mathbf{y}-\mathbf{X}\boldsymbol{\beta} \right\rVert^2 +\frac{1}{2\tau^2}
\left\lVert \boldsymbol{\beta} \right\rVert^2).
$$
したがって,事後分布のモード(MAP 推定量)は
\hat{\boldsymbol{\beta}}_{\mathrm{MAP}}= \arg\max_{\boldsymbol{\beta}}
\log p(\boldsymbol{\beta}\mid\mathbf{y},\sigma^2,\tau^2),
すなわち負の対数事後分布を最小化する問題として
\hat{\boldsymbol{\beta}}_{\mathrm{MAP}} = \arg\min_{\boldsymbol{\beta}}
\left\{
\left\lVert \mathbf{y} - \mathbf{X}\boldsymbol{\beta} \right\rVert^2 + \lambda \left\lVert \boldsymbol{\beta} \right\rVert^2
\right\},
\qquad \lambda = \frac{\sigma^2}{\tau^2},
と書くことができます。これはリッジ回帰の目的関数と一致しており,
リッジ推定量が「正規事前を仮定したベイズ回帰の MAP 推定」に
対応していることが分かります。
ちなみに,この目的関数は先ほど同様 $\boldsymbol{\beta}$ で微分することにより,
通常の閉形式解を得ることができます。目的関数を
$$
J(\boldsymbol{\beta})=
\left\lVert \mathbf{y} - \mathbf{X}\boldsymbol{\beta} \right\rVert^2 + \lambda \left\lVert \boldsymbol{\beta} \right\rVert^2
$$
とすれば、あとは
$$
\frac{\partial J(\boldsymbol{\beta})}{\partial \boldsymbol{\beta}} = -2\mathbf{X}^{\mathsf{T}}\mathbf{y} + 2\mathbf{X}^{\mathsf{T}}\mathbf{X}\boldsymbol{\beta} + 2\lambda \boldsymbol{\beta}.
$$
という感じで、
$$
\hat{\boldsymbol{\beta}} =
(\mathbf{X}^{\mathsf{T}}\mathbf{X} + \lambda \mathbf{I}_p)^{-1}
\mathbf{X}^{\mathsf{T}}\mathbf{y}
$$
がリッジ推定量(かつMAP推定量)となります。
結論
線形回帰において、
$\boldsymbol{\beta}$を定数とする最尤推定(MLE)はOLSと同義。
$\boldsymbol{\beta}$を正規事前とするMAP推定はRidge回帰と同義。
補足:ここでは $\boldsymbol{\beta}$ の分散を単一のスカラー $\tau^2$ としたが,
変数ごとに異なる分散 $\tau_j^2$ を許すこともできる。
その場合,階層モデルにより SNP ごとの効果分散を推定することになり,
BayesA/B/C 型の手法に一般化される。
参考文献
PRML: Pattern Recognition and Machine Learning
https://www.microsoft.com/en-us/research/wp-content/uploads/2006/01/Bishop-Pattern-Recognition-and-Machine-Learning-2006.pdf
BGLR: Genome-Wide Regression and Prediction with the BGLR Statistical Package
https://academic.oup.com/genetics/article/198/2/483/5935898