最尤推定法に基づく正規方程式の導出

被説明変数\(\boldsymbol{y}\)を、説明変数\(\boldsymbol{X}\)と回帰係数\(\boldsymbol{\beta}\)の線形結合と誤差項\(\boldsymbol{\varepsilon}\)によって表現する線形回帰モデル

\[ \newcommand{\b}[1]{\boldsymbol{#1}} \b{y} = \b{X} \b{\beta} + \b{\varepsilon} \]

を考える。また、以下を仮定する（古典的正規回帰モデル（Classical Normal Regression Model: CNRM） の仮定）

1. 説明変数\(\b{X}\)は非確率的である

2. \(E[Y] = \b{X \beta}\)であり、したがって誤差項の期待値はゼロ：\(E(\b{\varepsilon}) = 0\)

3. 誤差項の分散\(\sigma^2\)は一定（均一分散）であり、共分散はゼロ（独立性）：\(Var(\b{\varepsilon}) = \sigma^2 \b{I}\)

4. \(\b{X}\)の階数は\(k\)：\(rank(\b{X}) = k\)（\(\b{X^\top X}\)に逆行列が存在することの仮定）

5. \(Y\)は正規分布に従う（正規性）：\(Y \sim \mathcal{N}(\b{X\beta}, \sigma^2\b{I}), \quad \b{\varepsilon} \sim \mathcal{N}(0, \sigma^2\b{I})\)

特に3（標本がi.i.d.）と5（正規性）は最尤推定のために必要で、誤差項\(\varepsilon_i = y_i - \boldsymbol{x}_i^\top \boldsymbol{\beta}\)が従う分布型を仮定することで最尤推定が可能になる。

尤度としては、平均\(\boldsymbol{x}_i^\top \boldsymbol{\beta}\)の正規分布で観測値\(y_i\)が得られる確率\(\mathcal{N}(y_i| \boldsymbol{x}_i^\top \boldsymbol{\beta}, \sigma^2I)\)を用いて

\[ \newcommand{\b}[1]{\boldsymbol{#1}} L(\b{y}| \b{X}, \b{\beta}, \sigma) = \prod^N_{i=1} \mathcal{N}(y_i| \b{x}_i^\top \b{\beta}, \sigma^2I) \]

ただし、\(\mathcal{N}\)は正規分布

\[ \mathcal{N}(x|\mu, \sigma^2) = \frac{1}{\sqrt{2\pi}\sigma} \exp \left\{ - \frac{(x-\mu)^2}{2\sigma^2} \right\} , \hspace{1em} -\infty < x < \infty \]

対数尤度は

\[\begin{split} \begin{align} \ln L(\b{y}| \b{X}, \b{\beta}, \sigma) &= \sum^N_{i=1} \ln \mathcal{N}(y_i| \b{x}_i^\top \b{\beta}, \sigma^2) \\ &= - \frac{N}{2} \ln (2\pi) - \frac{N}{2} \ln \sigma^2 - \frac{1}{2\sigma^2} \sum^N_{i=1} (y_i - \b{x}_i^\top \b{\beta})^2 \end{align} \end{split}\]

導出メモ

まず尤度を軽く整理する

\[\begin{split} \begin{align} L(\mu, \sigma^2) &= \prod_{i=1}^n \frac{1}{\sqrt{2\pi\sigma^2}} \exp \left \{ - \frac{(x_i - \mu)^2}{2\sigma^2} \right \} \\ &= \left( \frac{1}{\sqrt{2\pi\sigma^2}} \right)^n \prod_{i=1}^n \exp \left \{ - \frac{(x_i - \mu)^2}{2\sigma^2} \right \} \end{align} \end{split}\]

対数をとると

\[\begin{split} \begin{align} \ln L(\mu, \sigma^2) &= n \ln \left( \frac{1}{\sqrt{2\pi\sigma^2}} \right) + \sum^n_{i=1} \ln \left( \exp \left\{ - \frac{(x_i - \mu)^2}{2\sigma^2} \right\} \right) \\ &= n \ln ( 2\pi\sigma^2 )^{-1/2} + \sum^n_{i=1} \left( - \frac{(x_i - \mu)^2}{2\sigma^2} \ln e \right) \\ &= - \frac{n}{2} \ln ( 2\pi\sigma^2 ) - \sum^n_{i=1} \frac{(x_i - \mu)^2}{2\sigma^2} \\ &= - \frac{n}{2} \ln (2\pi) - \frac{n}{2} \ln \sigma^2 - \frac{1}{2\sigma^2} \sum^n_{i=1} (x_i - \mu)^2 \end{align} \end{split}\]

（参考）対数の公式

\[\begin{split} \begin{align} \log 1 &= 0\\ \log_e e &= 1\\ \log x^a &= a\log x\\ \log \sqrt{a} &= \log a^{1/2} = \textstyle \frac{1}{2} \log a\\ \log ab &= \log a + \log b\\ \log \frac{a}{b} &= \log a - \log b \end{align} \end{split}\]

回帰係数の推定

対数尤度を\(\b{\beta}\)について微分した勾配を0と置いて解くと

\[\begin{split} \begin{align} \nabla \ln L(\b{y}| \b{X}, \b{\beta}, \sigma) &= \frac{1}{\sigma^2} \sum^N_{i=1} (y_i - \b{x}_i^\top \b{\beta})\b{x}_i^\top\\ &= \frac{1}{\sigma^2} (\sum^N_{i=1} y_i \b{x}_i^\top) - \b{\beta}^\top \frac{1}{\sigma^2} (\sum^N_{i=1} \b{x}_i \b{x}_i^\top)\\ &= \frac{1}{\sigma^2} (\b{X}^\top \b{y})^\top - \b{\beta}^\top \frac{1}{\sigma^2} \b{X}^\top \b{X}\\ &= 0\\ \to \hat{\b{\beta}} &= (\b{X}^\top \b{X})^{-1} \b{X}^\top \b{y}\\ \end{align} \end{split}\]

となり、最尤推定量と最小二乗推定量が同じ方程式になることがわかる

導出メモ

対数尤度を\(\b{\beta}\)について微分した勾配を整理すると

\[\begin{split} \begin{align} \nabla \ln L(\b{y}| \b{X}, \b{\beta}, \sigma) &= \frac{1}{\sigma^2} \sum^N_{i=1} (y_i - \b{x}_i^\top \b{\beta})\b{x}_i^\top\\ &= \frac{1}{\sigma^2} \sum^N_{i=1} y_i \b{x}_i^\top - \frac{1}{\sigma^2} \sum^N_{i=1} \b{x}_i^\top \b{\beta} \b{x}_i^\top\\ &= \frac{1}{\sigma^2} \sum^N_{i=1} y_i \b{x}_i^\top - \frac{1}{\sigma^2} \sum^N_{i=1} \b{\beta}^\top \b{x}_i \b{x}_i^\top\\ &= \frac{1}{\sigma^2} (\sum^N_{i=1} y_i \b{x}_i^\top) - \b{\beta}^\top \frac{1}{\sigma^2} (\sum^N_{i=1} \b{x}_i \b{x}_i^\top)\\ &= \frac{1}{\sigma^2} (\b{X}^\top \b{y})^\top - \b{\beta}^\top \frac{1}{\sigma^2} \b{X}^\top \b{X} \end{align} \end{split}\]

（\(\sum^N_{i=1} y_i \b{x}_i^\top = (\b{X}^\top \b{y})^\top\)は行ベクトルで\(\b{\beta}^\top \b{X}^\top \b{X}\)も行ベクトルなので演算できる）

これを0とおいて

\[ 0 = \frac{1}{\sigma^2} (\b{X}^\top \b{y})^\top - \b{\beta}^\top \frac{1}{\sigma^2} \b{X}^\top \b{X} \]

として整理すると

\[ \b{\beta}^\top \frac{1}{\sigma^2} \b{X}^\top \b{X} = \frac{1}{\sigma^2} (\b{X}^\top \b{y})^\top \]

\[\begin{split} \begin{align} &\to \b{\beta}^\top \frac{1}{\sigma^2} \b{X}^\top \b{X} = \frac{1}{\sigma^2} (\b{X}^\top \b{y})^\top\\ &\to \b{\beta}^\top (\b{X}^\top \b{X}) = (\b{X}^\top \b{y})^\top\\ &\to \b{\beta}^\top = (\b{X}^\top \b{X})^{-1} (\b{X}^\top \b{y})^\top\\ &\to \b{\beta} = (\b{X}^\top \b{X})^{-1} \b{X}^\top \b{y}\\ \end{align} \end{split}\]

（\((\b{X}^\top \b{X})\)は対角行列なので\((\b{X}^\top \b{X})^\top = (\b{X}^\top \b{X})\)）

※内積の結果はスカラーになるので、\(\b{x}^\top \b{\beta} = \b{\beta}^\top \b{x}\)である

2つの行列\(A, B\)があるとき、その行ベクトル\(a_i, b_i\)の外積\(\boldsymbol{a}_i \boldsymbol{b}_i^\top\)と行列の積は一致する

\[ A^\top B = \sum^n_{i=1} \boldsymbol{a}_i \boldsymbol{b}_i^\top \]

ベクトル\(\boldsymbol{b} = (b_1, b_2, \dots, b_n)^\top\)とベクトル\(\boldsymbol{a}_i = (a_{i1}, a_{i2}, \cdots, a_{im})\)があるとき、スカラーの\(b_i\)との積の和は

\[\begin{split} \begin{align} &\sum_{i=1}^n b_i \boldsymbol{a}_i^\top\\ &= \sum_{i=1}^n (b_i a_{i1}, \cdots, b_i a_{im})\\ &= \begin{pmatrix} b_1 a_{11} & \cdots & b_1 a_{1m}\\ + & & +\\ \vdots & \cdots & \vdots\\ + & & +\\ b_n a_{n1} & \cdots & b_n a_{nm}\\ \end{pmatrix}\\ \end{align} \end{split}\]

ベクトル\(\b{a}_i\)を集めた行列\(A\)を考えたとき、

\[\begin{split} \b{A}^\top \b{b} = \begin{pmatrix} b_1 a_{11} + \cdots + b_n a_{n1}\\ \vdots\\ b_1 a_{m1} + \cdots + b_n a_{nm}\\ \end{pmatrix} \end{split}\]

ゆえに

\[ \sum_{i=1}^n b_i \boldsymbol{a}_i^\top = (A^\top b)^\top \]

導出メモ

スカラーの場合、合成関数の微分より

\[ \frac{d (y - \hat{y})^2 }{d \hat{y}} = 2 (y - \hat{y}) \times (-1) \]

の形になるがベクトルなので

\[ \frac{\partial (y_i - \b{x}_i^\top \b{\beta})^2 }{\partial \beta_j} = - 2 (y_i - \b{x}_i^\top \b{\beta}) \b{x}_i^\top \]

これを各\(j\)について計算してベクトル（勾配）とし、それを\(i\)ごとに和する

標準偏差\(\sigma\)の導出

\[\begin{split} \begin{align} \frac{\partial \ln L(\b{y}| \b{X}, \b{\beta}, \sigma) }{ \partial \sigma } &= \frac{\partial}{ \partial \sigma } \left(- \frac{N}{2} \ln (2\pi) - \frac{N}{2} \ln \sigma^2 - \frac{1}{2\sigma^2} \sum^N_{i=1} (y_i - \b{x}_i^\top \b{\beta})^2 \right)\\ &= - \frac{N}{2} \frac{1}{\sigma^2} 2\sigma + \frac{4\sigma}{4\sigma^4} \sum^N_{i=1} (y_i - \b{x}_i^\top \b{\beta})^2\\ &= - \frac{N}{\sigma} + \frac{1}{\sigma^3} \sum^N_{i=1} (y_i - \b{x}_i^\top \b{\beta})^2\\ \end{align} \end{split}\]

これをゼロにする\(\sigma^2\)は

\[\begin{split} \begin{align} &- \frac{N}{\sigma} + \frac{1}{\sigma^3} \sum^N_{i=1} (y_i - \b{x}_i^\top \b{\beta})^2 = 0\\ &\to \frac{N}{\sigma} = \frac{1}{\sigma^3} \sum^N_{i=1} (y_i - \b{x}_i^\top \b{\beta})^2\\ &\to \frac{\sigma^3}{\sigma} = \frac{\sigma^3}{N \sigma^3} \sum^N_{i=1} (y_i - \b{x}_i^\top \b{\beta})^2\\ &\to \sigma^2 = \frac{1}{N} \sum^N_{i=1} (y_i - \b{x}_i^\top \b{\beta})^2\\ \end{align} \end{split}\]

となる

\(\b{\beta}\)の最尤推定量 \(\hat{\b{\beta}}_{ML} = (\b{X}^\top \b{X})^{-1} \b{X}^\top \b{y}\)を代入すると

\[\begin{split} \begin{align} \hat{\sigma}^2_{ML} &= \frac{1}{N} \sum^N_{i=1} (y_i - \b{x}_i^\top \hat{\b{\beta}}_{ML})^2\\ &= \frac{ \b{\varepsilon}^\top \b{\varepsilon} }{N}\\ \end{align} \end{split}\]