IRLS で解く最尤推定量

一般化線形モデルの対数尤度関数から最尤推定量を計算する際に \(\frac{\partial l(\beta)}{\partial \beta}=0\) の解を求める必要がある。この方程式をニュートン・ラフソン法により求めることができる。

ニュートン・ラフソン法

関数 y = f(x) があり f(x)= 0 となる解を近似法によって求めるのがニュートン・ラフソン法である。

ニュートンラフソン法

f(x) = 0 の解を求めるために、グラフ上の任意の点 x = x1 における接線(緑色)を計算する。次に、x = x1 における接線と x 軸の交点を x2 として、x = x2 における接線(青色)を計算する。さらに、次に、x = x2 における接線と x 軸の交点を x3 として、x = x3 における接線(赤色)を計算する。このように x4、x5 …と繰り返して求めていけば、xn が限りなく f(x) = 0 の解に近づく。このようなアルゴリズムを用いた方程式の解法をニュートン・ラフソン法という。

xm から xm+1 を求めるとき、以下のように数式化することができる。

\[ 0 - f(x_{m}) = f'(x_{m})(x_{m+1} - x_{m}) \] \[ \therefore 0 = f(x_{m}) + f'(x_{m})(x_{m+1} - x_{m}) \]

ニュートン・ラフソン法(IRLS)による最尤推定量の求め方

ニュートン・ラフソン法を応用するために、上で説明した関数 f を \( \frac{\partial l(\beta)}{\partial \beta} \) とみなせばよい。(m+1) 回目のパラメータ β(m+1) を計算する際に以下のように計算する。

\[ 0 = \left . \frac{\partial l}{\partial \beta} \right \vert_{(m)} + \left . \frac{\partial ^{2} l}{\partial \beta \partial \beta^{T}} \right \vert_{(m)} (\beta^{(m+1)}-\beta^{(m)}) \] \[ \therefore \beta^{(m+1)} = \beta^{(m)} - \frac{\left . \frac{\partial l}{\partial \beta} \right \vert_{(m)}}{\left . \frac{\partial ^{2} l}{\partial \beta \partial \beta^{T}}\right \vert_{(m)} } \]

2 階微分の計算が難しいため、その期待値を利用する。すなわち、\( \frac{\partial ^{2} l}{\partial \beta \partial \beta^{T}} \) の代わりに \( E\left[ \frac{\partial ^{2} l}{\partial \beta \partial \beta^{T}} \right]\) を用いる。

このように対数尤度関数の 2 階微分の期待値を利用した方法は、反復重み付け最小二乗法(IRLS)とも呼ばれている。

\[ \begin{eqnarray} \beta^{(m+1)} &=& \beta^{(m)} - \frac{\left . \frac{\partial l}{\partial \beta} \right \vert_{(m)}}{\left .E\left[ \frac{\partial ^{2} l}{\partial \beta \partial \beta^{T}} \right]\right \vert_{(m)} } \\ &=& \beta^{(m)} + \frac{\left . \frac{\partial l}{\partial \beta} \right \vert_{(m)}}{\left . - E\left[ \frac{\partial ^{2} l}{\partial \beta \partial \beta^{T}} \right]\right \vert_{(m)} } \\ &=& \beta^{(m)} + \frac{\mathbf{U}^{m}}{\mathbf{I}^{(m)}} \\ &=& \beta^{(m)} + \left [ \mathbf{I}^{(m)} \right ] ^{-1} \mathbf{U}^{m} \\ &=& \beta^{(m)} + \left [ \mathbf{X}^{T}\mathbf{W}\mathbf{X}^{(m)} \right ] ^{-1} \mathbf{U}^{m} \\ \end{eqnarray} \]

この式を利用して、ニュートン・ラフソン法により、\( l(\beta^{(n+1)}) - l(\beta^{(n)}) \lt \delta \) になるまで繰り返す(\(\delta \rightarrow 0\))。しかし、観測されたデータの分布に偏りがある場合は、計算が収束しない場合がある。例えば、ゼロ過剰の二項分布などがそうである。

References

  1. P.M.E.Altham, Statistical Laboratory, University of Cambridge. Introduction to Generalized Linear Modelling. 2011. PDF
  2. Dobson AJ. An Introduction to Generalized Linear Models. Second Edition. 2002.
  3. ニュートン法. Wikipedia