生物統計学

対数尤度関数の最尤推定

確率変数 Y を一般化線形モデルでモデル化するとき、パラメーターを β とし、デザイン行列を X とし、リンク関数を g とすると、モデルは次式のようにかける。このとき、一般化線形モデルに組み込まれたパラメーター β を最尤法により推定する手順を示す。

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

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

確率変数 Y の密度関数が、次式のような指数型分布族で表せるものとする。

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

密度関数が f(Y; θ) のとき、その尤度関数 L(θ Y) = f(Y; θ) も同様な式で表すことができ、このとき対数尤度関数は以下のようにかける。

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

いま、n 個の確率変数 y1, ..., yn があり、このときの対数尤度関数は以下のように書ける。最適なパラメーター β を求めるには、この対数尤度関数を最大にする β を求めればよい。

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

対数尤度関数の 2 階微分は負の値を取る（証明）ため、対数尤度関数が最大値をとるときは、1 階微分が 0 になるときである。

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

対数尤度関数の 1 階微分が 0 となる β を求めるのに、ニュートン・ラフソン法を利用する。ニュートン・ラフソン法を利用するには、1 階微分だけでなく、2 階微分微分も必要である。ここで、まず対数尤度関数の 1 階微分と 2 階微分を求めていく。

最尤推定量（対数尤度関数の 1 階微分）

対数尤度関数の 1 階微分は次のように書ける。

\[\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 \mu_{i}}{\partial \theta_{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 項については、リンク関数の両辺を β で微分すればよい。ここで、xi· をデザイン行列 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} \]

以上により、対数尤度関数を β で微分した式は、以下のように書き換えることができる。対数尤度関数を微分した量をスコア関数はスコア関数として呼ばれることもある。スコア関数は一般的に記号 U で表すことが多い。

\[ \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} \]

対数尤度関数の 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] \]

右辺は対数尤度関数を β で微分したものであるので、対数尤度関数の 1 階微分で得られた結果を代入して整理する。

\[ \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 階導関数をスコア関数といい、スコア関数の分散をフィッシャー情報量という。パラメーターが p 個あるとき、p 個のパラメーターそれぞれに対するフィッシャー情報量が求まる。これらのフィッシャー情報量を p×p の正定値対称行列で表した行列をフィッシャー情報行列という。フィッシャー情報量行列は 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 アルゴリズム）を利用して、対数尤度関数を最大にする β を求められるようになる。

ニュートン・ラフソン法

求めたいものは \( \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)}} \]

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.