尤度方程式
尤度方程式(ゆうどほうていしき、likelihood equation)とは、統計学において、対数尤度関数の極値条件を与える方程式の事。統計的推定法の一つである最尤法において、尤度関数を最大化する最尤推定値を求める際に用いられる。
X1, …, Xnを同一の確率分布に従う独立な確率変数とし、x(x1,…, xn)Tをその観測値とする。ここで確率分布は確率密度関数f(x! θ)をもつ連続分布とする。但し、θ(θ1, .., θp)Tはp個の分布のパラメータを表すものとし、パラメータ空間Θ ⊂ Rpに値を持つものとする。このとき、尤度関数L(θ,x)と対数尤度関数l(θ,x)は以下で定義される。
L(\boldsymbol{\theta}, \mathbf{x})=\prod_{i=1}^n f(x_i|\boldsymbol{\theta})
l(\boldsymbol{\theta}, \mathbf{x})=\ln{L(\boldsymbol{\theta}, \mathbf{x})}=\sum_{i=1}^n \ln{f(x_i|\boldsymbol{\theta})}
最尤法では、対数尤度関数l(θ,x)を最大化するθが最尤推定値として定まる。このとき、最尤推定値θは、l(θ,x)に対する次の極値条件を満たす。
\frac{\partial}{\partial \boldsymbol{\theta}l(\boldsymbol{\theta}, \mathbf{x})=\mathbf{0}
θに関するこの方程式を尤度方程式という。左辺の勾配ベクトル
\mathbf{S}(\mathbf{x},\boldsymbol{\theta}):=\frac{\partial}{\partial \boldsymbol{\theta}l(\boldsymbol{\theta}, \mathbf{x})
は、スコア関数、もしくは単にスコアと呼ばれる。多くの場合、最尤推定値の推定は、尤度方程式を解く問題、すなわち、スコアをゼロとするパラメータθ∈ Θを求める問題に帰着する。
Xi (i1,..,n)が平均をμ、分散をσ2とする正規分布に従うとする(X ∼ N(μ, σ2))。このとき、対数尤度関数は
l(\mu,\sigma^2,\mathbf{x})=-\frac{n}{2} \ln{2\pi}-\frac{n}{2} \ln{\sigma^2}-\frac{1}{2\sigma^2}\sum_{i=1}^n(x_i-\mu)^2
であり、尤度方程式は
\frac{\partial l(\mu,\sigma^2,\mathbf{x})}{\partial \mu}=\frac{1}{\sigma^2}\sum_{i=1}^n(x_i-\mu)=0
\frac{\partial l(\mu,\sigma^2,\mathbf{x})}{\partial \sigma^2}=
-\frac{n }{2\sigma^2}+\frac{1}{2(\sigma^2)^2}\sum_{i=1}^n(x_i-\mu)^2=0
となる。これらを整理すると最尤推定値として
\hat{\mu}=\frac{1}{n} \sum_{i=1}^nx_i
\hat{\sigma^2}=\frac{1}{n}\sum_{i=1}^n(x_i-\mu)^2
を得る。
Xi (i1,..,n)が形状パラメータをβ、尺度パラメータをηとするワイブル分布に従うとする。このとき、対数尤度関数は
l(\eta,\beta,\mathbf{x})=n \ln{\beta}-n \beta \ln{\eta} +(\beta-1)\sum_{i=1}^n \ln{x_i}-\frac{1}{\eta^\beta}\sum_{i=1}^nx_i^{\beta}
であり、尤度方程式は
\frac{\partial l(\eta,\beta,\mathbf{x})}{\partial \eta}=-\frac{n\beta}{\eta}- \frac{\beta}{\eta^{(\beta+1)}\sum_{i=1}^nx_i^{\beta}=0
\frac{\partial l(\eta,\beta,\mathbf{x})}{\partial \beta}=\frac{n}{\beta}-n\ln{\eta}+\sum_{i=1}^n\ln{x_i}+ \frac{\ln{\eta}{\eta^{\beta}\sum_{i=1}^nx_i^{\beta}
+ \frac{1}{\eta^{\beta}\sum_{i=1}^n\ln{x_i}x_i^{\beta}=0
となる。
これらを整理すると最尤推定値η、βが満たすべき関係式
\hat{\eta}=\left( \frac{1}{n}\sum_{i=1}^n x_i^{\hat{\beta} \right)^{\frac{1}{\hat{\beta}}
\frac{1}{\hat{\beta}+
\frac{1}{n}\sum_{i=1}^n\ln{x_i}-\frac{\sum_{i=1}^nx_i^{\hat{\beta}\ln{x_i} }{\sum_{i=1}^nx_i^{\hat{\beta}}=0
を得る。第二式を満たすβを数値的に求めれば、第一式よりηも定まる。
Xi (i1,..,n)が形状パラメータをα、尺度パラメータをβとするガンマ分布に従うとする(X ∼ Γ(α, β))。このとき、対数尤度関数は
l(\alpha,\beta,\mathbf{x})=-n \ln{\Gamma(\alpha)}-n \alpha \ln{\beta} +(\alpha-1)\sum_{i=1}^n \ln{x_i}-\frac{1}{\beta}\sum_{i=1}^nx_i
であり、尤度方程式は
\frac{\partial l(\alpha,\beta,\mathbf{x})}{\partial \alpha}=-n \psi(\alpha)-n \ln{\beta} +(\alpha-1)\sum_{i=1}^n \ln{x_i}=0
\frac{\partial l(\alpha,\beta,\mathbf{x})}{\partial \beta}=
-\frac{n \alpha}{\beta}+\frac{1}{\beta^2}\sum_{i=1}^nx_i=0
となる。
ここではψ(α)はガンマ関数の対数微分であるディガンマ関数を表す。これらを整理すると最尤推定値β、αが満たすべき関係式
\hat{\beta}=\frac{1}\hat\alpha \frac{1}{n}\sum_{i=1}^nx_i
\hat{\alpha}=\frac{\frac{1}{n}\sum_{i=1}^nx_i}{\left (\prod_{i=1}^nx_i \right)^{\frac{1}{n} }\exp{(\psi(\hat{\alpha}))}
を得る。第二式を満たすαを数値的に求めれば、第一式よりβも定まる。
尤度方程式が解析的に解けない場合、S(θ*)0を満たすθ*∈ Θを数値的に求めることが必要となる。
ニュートン=ラフソン法では、反復計算により、最適解θ*を求める。反復計算のkステップ目で求まったパラメータをθ(k)とする。スコア関数はテイラー展開により、
\mathbf{S}(\mathbf{x},\boldsymbol{\theta})\simeq
\mathbf{S}(\mathbf{x},\boldsymbol{\theta}^(k))-I(\boldsymbol{\theta}^{(k)})(\boldsymbol{\theta}-\boldsymbol{\theta}^{(k)})
と一次近似できる。ここでI(θ)は、
I(\boldsymbol{\theta})=
-\frac{\partial^2}{\partial \boldsymbol{\theta}\partial \boldsymbol{\theta}^T}\ln{L(\boldsymbol{\theta}, \mathbf{x})}
で与えられる、対数尤度関数のヘッセ行列の符号を変えた行列である。ニュートン=ラフソン法では、左辺をゼロとおくことで、θ(k+1)を与える更新式
\boldsymbol{\theta}^{(k+1)}
=\boldsymbol{\theta}^{(k)}+
I(\boldsymbol{\theta}^{(k)})^{-1}\mathbf{S}(\mathbf{x},\boldsymbol{\theta}^(k))
を定める。
ニュートン=ラフソン法は、最適解θ*の近傍で二次収束するため、収束が早い。すなわち、θ*の十分近くの適切な初期値を与えれば、
||\boldsymbol{\theta}^(k)-\boldsymbol{\theta}^{\ast}||
\leq K ||\boldsymbol{\theta}^(k)-\boldsymbol{\theta}^{\ast}||^2
を満たす正の定数Kが存在する。
一方で、ニュートン=ラフソン法は各ステップで、対数尤度関数のヘッセ行列から定まるI(θ)の逆行列を計算する、もしくは、p次の連立方程式を解くことが必要となる。これらの計算量はO(p3)のオーダーであり、パラメータ数pが増えると、計算負荷が急激に増える。また、初期値の設定によっては、I(θ)は正定値とはならず、最適解θ*に収束しない場合がある。
ニュートン=ラフソン法においては、各ステップで負の対数尤度関数の二階微分であるI(θ)を計算する必要がある。このI(θ)を求める計算は、場合によっては煩雑となる。分布によっては、I(θ)の期待値であるフィッシャー情報行列
J(\boldsymbol{\theta})=
E_{\boldsymbol{\theta} \left [-\frac{\partial^2}{\partial \boldsymbol{\theta}\partial \boldsymbol{\theta}^T}\ln{L(\boldsymbol{\theta}, \mathbf{x})} \right ]
=E_{\boldsymbol{\theta} \left [ \frac{\partial}{\partial \boldsymbol{\theta} \ln{L(\boldsymbol{\theta}, \mathbf{x}) \frac{\partial}{\partial \boldsymbol{\theta}^T} \ln{L(\boldsymbol{\theta}, \mathbf{x})} \right ]
が、より簡潔に求まるため、I(θ)をJ(θ)で代用し、反復計算を
\boldsymbol{\theta}^{(k+1)}
=\boldsymbol{\theta}^{(k)}+
J(\boldsymbol{\theta}^{(k)})^{-1}\mathbf{S}(\mathbf{x},\boldsymbol{\theta}^(k))
とする。この方法をフィッシャーのスコア法と呼ぶ。
フィッシャー情報行列は非負定値であるため、ニュートン=ラフソン法でのI(θ)の正定値性の問題を回避することができる。