対数尤度関数をニュートン・ラフソン法で解く

対数尤度関数の最尤推定

一般化線形モデルの対数尤度関数

一般化線形モデルに組み込まれたパラメーターを β を最尤法により推定する手順を示す。まず、指数型分布族の確率質量・密度関数を以下の形で表すものとする。モデルを立てる際にデザイン行列を X で表すので、ここでは確率変数を Y とする。

\[ f(Y; \theta) = h(Y)\exp \left( \eta (\theta) Y - A(\theta) \right) \]

その対数尤度関数は以下のように書ける。

\[ l(\theta; Y) = \log f(Y; \theta) = \log h(Y) + \eta (\theta) Y - A(\theta) \]

ここで、n 個の確率変数があり、それらに対する対数尤度関数は以下のように書ける。

\[ l(\mathbf{\theta}; \mathbf{y}) = \sum \log h(y_{i}) + \sum \eta (\theta_{i}) y_{i} - \sum A(\theta_{i}) \]

また、一般化線形モデルを組み立てる際にデザイン行列を X、リンク関数を g とする。すなわち、

\[ E[y_{i}] = \mu_{i} \] \[ g(\mu_{i}) = \mathbf{x}_{i}^{T}\mathbf{\beta} \]

最尤推定量(対数尤度関数の 1 階微分)

対数尤度関数の 2 階微分は負の値を取る(証明)ため、対数尤度関数が最大値をとるときは、すなわち

\[\frac{\partial l}{\partial \beta} = 0\]

のときである。まず、左辺の対数尤度関数の微分を計算しやすいように式変換する。

微分の連鎖率により、\(\frac{\partial l}{\partial \beta}\) は下のように 3 つの項に分解することができる。

\[ \frac{\partial l}{\partial \beta} = \sum _{i=1}^{n} \frac{\partial l}{\partial \theta_{i}} \frac{\partial \theta_{i}}{\partial \mu_{i}}\frac{\partial \mu_{i}}{\partial \beta} \]

第 1 項については対数尤度関数を θi で微分すれば解が得られる。

\[ \begin{eqnarray} \frac{\partial l}{\partial \theta_{i}} &=& \eta'(\theta_{i})y_{i} - A'(\theta_{i}) \\ &=& \eta'(\theta_{i})\left(y_{i} - \frac{A'(\theta_{i})}{\eta'(\theta_{i})}\right) \\ &=& \eta'(\theta_{i})\left(y_{i} - \mu_{i}\right) \end{eqnarray} \]

第 2 項については、\( \frac{A'(\theta)}{\eta'(\theta)}=\mu \) の関係を利用すれば解が得られる。

\[ \frac{\partial \theta_{i}}{\partial \mu_{i}} = \frac{A''(\theta_{i})\eta'(\theta_{i}) - A'(\theta_{i})\eta''(\theta_{i})}{(\eta'(\theta_{i}))^{2}} = \eta'(\theta_{i})Var[y_{i}] \]

第 3 項については、リンク関数の両辺を β で微分すればよい。ここで、x をデザイン行列 X の i 行目の行ベクトルとする。

\[ \begin{eqnarray} \frac{\partial}{\partial \mathbf{\beta}}g(\mu_{i}) &=& \frac{\partial}{\partial \mathbf{\beta}} \mathbf{x}_{i\cdot}\mathbf{\beta} \\ \frac{\partial \mu_{i}}{\partial \mathbf{\beta}}\frac{\partial }{\partial \mu_{i}} g(\mu_{i}) &=& \mathbf{x}_{i\cdot} \\ \therefore \frac{\partial \mu_{i}}{\partial \mathbf{\beta}} &=& \frac{\mathbf{x}_{i\cdot}^{T} }{g'(\mu_{i})} \end{eqnarray} \]

そこで、対数尤度関数を β で微分した数式は以下のように書き換えることができる。

\[ \begin{eqnarray} \frac{\partial l}{\partial \beta} &=& \sum _{i=1}^{n} \left( \left(\eta'(\theta_{i})\left(y_{i} - \mu_{i} \right) \right) \left( \frac{1}{\eta'(\theta_{i})Var[y_{i}]}\right) \left(\frac{\mathbf{x}_{i}^{T} }{g'(\mu_{i})} \right)\right)\\ &=& \sum _{i=1}^{n} \frac{(y_{i} - \mu_{i})\mathbf{x}^{T}_{i \cdot}}{g'(\mu_{i})Var[y_{i}]} = \mathbf{U} \end{eqnarray} \]

対数尤度関数を微分した量をスコア関数ともよばれている。スコア関数は U として記述されることがある。

これで、対数尤度関数の 1 階微分を計算しやすい形に式変換を行った。次にこの関数が 0 を取るときの解を計算する。これにはニュートン・ラフソン法がより利用される。

ニュートン・ラフソン法

求めたいものは \( \frac{\partial l}{\partial \beta} = U = 0 \) である。そこで、対数尤度関数の 1 階微分を第 2 項までテイラー展開行う。

\[ \mathbf{0} \simeq \frac{\partial l}{\partial \mathbf{\beta} ^{(0)}} + \frac{\partial ^{2} l}{\partial \mathbf{\beta}^{(0)} \partial \mathbf{\beta} ^{(0)^{T}} } (\mathbf{\beta} - \mathbf{\beta}^{(0)}) \]

これを整理すると以下のようになる。

\[ \mathbf{\beta} \simeq \mathbf{\beta}^{(0)} - \left \{ \frac{\partial ^{2} l}{\partial \mathbf{\beta}^{(0)} \partial \mathbf {\beta} ^{(0)^{T}} } \right \}^{-1} \frac{\partial l}{\partial \mathbf{\beta} ^{(0)}} \]

ニュートン・ラフソン法ではこの式を繰り返し用いて、最適な β を計算する(β を最適化する)。

\[ \mathbf{\beta}^{(m+1)} \simeq \mathbf{\beta}^{(m)} - \left \{ \frac{\partial ^{2} l}{\partial \mathbf{\beta}^{(m)} \partial \mathbf {\beta} ^{(m)^{T}} } \right \}^{-1} \frac{\partial l}{\partial \mathbf{\beta} ^{(m)}} \]

このようにニュートン・ラフソン法では、対数尤度関数の 2 階微分を利用するので、対数尤度関数の 2 階微分を求める。

対数尤度関数の 2 階微分

\( \frac{\partial ^{2}l}{\partial \beta \partial \beta ^{T}} \) について見ていく。

まず、分布関数については \( \int f(y; β) dy = 1 \) が成り立つので、両辺を β および βT で微分する。(以下の式変形において l(β y) = logf(y; β) を利用した部分がある)
\[ \begin{eqnarray} 0 &=& \frac{\partial ^{2}}{\partial \beta \partial \beta ^{T}} \int f(y; \beta)dy \\ &=& \frac{\partial}{\partial \beta} \left( \frac{\partial}{\partial \beta^{T}} \int \exp l(\beta; y)dy \right) \\ &=& \frac{\partial}{\partial \beta} \left( \int \left(\frac{\partial}{\partial \beta^{T}} l(\beta; y)\right) \exp l(\beta; y)dy \right) \\ &=& \int\frac{\partial}{\partial \beta} \frac{\partial}{\partial \beta^{T}} l(\beta; y) \exp l(\beta; y)dy \\ && + \int \left(\frac{\partial}{\partial \beta^{T}} l(\beta; y)\right) \left(\frac{\partial}{\partial \beta} l(\beta; y)\right) \exp l(\beta; y)dy \\ &=& \int\frac{\partial}{\partial \beta} \frac{\partial}{\partial \beta^{T}} l(\beta; y) f(y;\beta;)dy \\ && + \int \left(\frac{\partial}{\partial \beta^{T}} l(\beta; y)\right) \left(\frac{\partial}{\partial \beta} l(\beta; y)\right) f(y;\beta)dy \\ &=& E\left[ \frac{\partial ^{2}}{\partial \beta \partial \beta^{T}} l(\beta; y) \right] + E\left[ \frac{\partial}{\partial \beta^{T}} l(\beta; y) \frac{\partial}{\partial \beta} l(\beta; y) \right] \end{eqnarray} \]

よって、

\[ \therefore E\left[ \frac{\partial ^{2}}{\partial \beta \partial \beta^{T}} l(\beta; y) \right] = - E\left[ \frac{\partial}{\partial \beta} l(\beta; y) \frac{\partial}{\partial \beta ^{T}} l(\beta; y) \right] \]

右辺は対数尤度関数を β で微分したものであるので、上で計算したものを代入して整理する。

\[ \begin{eqnarray} E\left[ \frac{\partial ^{2}}{\partial \beta \partial \beta^{T}} l(\beta; y) \right] &=& -E\left[ \sum_{i=1}^{n} \frac{(y_{i} - \mu_{i})^{2} \mathbf{x}^{T}_{i\cdot}\mathbf{x}_{i\cdot}}{(g'(\mu_{i})Var[y_{i}])^{2}} \right] \\ &=& - \sum_{i=1}^{n} \frac{E\left[ (y_{i} - \mu_{i})^{2} \right]\mathbf{x}_{i\cdot}^{T}\mathbf{x}_{i\cdot} }{(g'(\mu_mu)Var[y_{i}])^{2}} \\ &=& - \sum_{i=1}^{n} \frac{Var[y_{i}] \mathbf{x}_{i\cdot}^{T}\mathbf{x}_{i\cdot} }{(g'(\mu_{i})Var[y_{i}])^{2}} \\ &=& - \sum_{i=1}^{n} \frac{\mathbf{x}_{i\cdot}^{T}\mathbf{x}_{i\cdot}}{(g'(\mu_{i}))^{2} Var[y_{i}]} \\ &=& - \sum_{i=1}^{n} w_{i} \mathbf{x}_{i\cdot}^{T}\mathbf{x}_{i\cdot} \end{eqnarray} \]

ただし、

\[ w_{i} = \frac{1}{(g'(\mu_{i}))^{2} Var[y_{i}]} \\ \]

さらに W を w からなる行列とし、すなわち、

\[ \mathbf{W} = \begin{pmatrix} w_{1} & 0 & 0 & \cdots & 0 \\ 0 & w_{2} & 0 & \cdots & 0 \\ \vdots & \vdots & \vdots & \ddots & \vdots \\ 0 & 0 & 0 & \cdots & w_{n} \\ \end{pmatrix} \]

このとき、以下の式が成り立つ。

\[ E\left[ \frac{\partial ^{2}}{\partial \beta \partial \beta^{T}} l(\beta; y) \right] = -\mathbf{X}^{T}\mathbf{W}\mathbf{X} \]

対数尤度関数の 1 階導関数はスコア関数と呼ばれる。スコア関数の分散をフィッシャー情報量という。ββT の共分散を表し、フィッシャー情報行列となっている(計算)。フィッシャー情報量行列は I または \(\mathfrak{F}\) (ドイツ文字の F)と記述されるため、上式は以下のように記述してもよい。

\[ I_{\mathbf{y}}(\beta) = \mathfrak{F} = - E\left[ \frac{\partial ^{2}}{\partial \beta \partial \beta^{T}} l(\beta; y) \right] = \mathbf{X}^{T}\mathbf{W}\mathbf{X} \]

ここで求めた XTWX が対数尤度関数の 2 階微分の期待値にあたる。これで 1 階微分と 2 階微分の期待値が求められたので、ここで始めてニュートン・ラフソン法(IRLS アルゴリズム)が利用できるようになる。

References

  1. James WH, Joseph MH. Generalized Linear Models and Extensions. Third Edition. 2012.
  2. P.M.E.Altham, Statistical Laboratory, University of Cambridge. Introduction to Generalized Linear Modelling. 2011. PDF
  3. Dobson AJ. An Introduction to Generalized Linear Models. Second Edition. 2002.