対数尤度関数の 1 階微分を計算しやすい形に式変換しニュートン・ラフソン法により解く

対数尤度関数の最尤推定量

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

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

最尤推定量

対数尤度関数の 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 階微分の値をそのまま用いるのではなく、その期待値を利用する場合もある。

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.