尤度方程式(ゆうどほうていしき、英: likelihood equation)とは、統計学において、対数尤度関数の極値条件を与える方程式の事。統計的推定法の一つである最尤法において、尤度関数を最大化する最尤推定値を求める際に用いられる。
概要
独立同分布を満たす
個の確率変数
とその観測値
を定義する。すなわち真の分布から
個の観測値(データ)が無作為抽出された状況を考える。
ここで確率密度関数
に従う確率モデルを導入する。ここで
は分布パラメータ群であり、パラメータ空間Θ ⊂ Rpに値を持つ。この確率モデルが
を最も良く説明する
を求めたい。ゆえに最尤推定をおこなう。
このとき独立同分布条件により、尤度関数
と対数尤度関数
は以下で定義される。
![{\displaystyle L({\boldsymbol {\theta }}|{\boldsymbol {d}})=\prod _{i=1}^{n}f(X=d_{i}|{\boldsymbol {\theta }})}](https://wikimedia.org/api/rest_v1/media/math/render/svg/b68a437a72841597eb4f704210d42d5c6ad1afde)
![{\displaystyle l({\boldsymbol {\theta }}|{\boldsymbol {d}})=\ln {L({\boldsymbol {\theta }}|{\boldsymbol {d}})}=\sum _{i=1}^{n}\ln {f(X=d_{i}|{\boldsymbol {\theta }})}}](https://wikimedia.org/api/rest_v1/media/math/render/svg/74dd39c7ea9890273a4f5f83035f96f7485b932e)
すなわちあるデータ群に対するモデルの尤度関数は、各観測値に対する尤度関数の積(対数尤度の場合は和)となる。
最尤法では対数尤度関数を最大化する
が最尤推定値
として定まる。このとき
は次の極値条件を満たす。
![{\displaystyle {\frac {\partial }{\partial {\boldsymbol {\theta }}}}l({\boldsymbol {\theta }}|{\boldsymbol {d}})=\mathbf {0} }](https://wikimedia.org/api/rest_v1/media/math/render/svg/41221008f751acc53eaf9267d60176003551c74a)
この方程式を尤度方程式という。左辺の勾配ベクトル:
![{\displaystyle \mathbf {S} ({\boldsymbol {d}},{\boldsymbol {\theta }}):={\frac {\partial }{\partial {\boldsymbol {\theta }}}}l({\boldsymbol {\theta }}|{\boldsymbol {d}})}](https://wikimedia.org/api/rest_v1/media/math/render/svg/3db615228b7a62d30af79cc61f8d3b95733230b3)
は、スコア関数、もしくは単にスコアと呼ばれる。多くの場合、最尤推定値の推定は、尤度方程式を解く問題、すなわち、スコアをゼロとするパラメータθ∈ Θを求める問題に帰着する。
例
正規分布
Xi (i=1,..,n)が平均をμ、分散をσ2とする正規分布に従うとする(X ∼ N(μ, σ2))。このとき、対数尤度関数は
![{\displaystyle 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}}](https://wikimedia.org/api/rest_v1/media/math/render/svg/3cbdf21d9ac6e18797a6b471c70a206b58085311)
であり、尤度方程式は
![{\displaystyle {\frac {\partial l(\mu ,\sigma ^{2},\mathbf {x} )}{\partial \mu }}={\frac {1}{\sigma ^{2}}}\sum _{i=1}^{n}(x_{i}-\mu )=0}](https://wikimedia.org/api/rest_v1/media/math/render/svg/d64e63684de7aa3c167dbf8644d9f156a9f1b88b)
![{\displaystyle {\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}](https://wikimedia.org/api/rest_v1/media/math/render/svg/70ac493f8baf6116b73ff8c592d9ae4aa94f66d8)
となる。これらを整理すると最尤推定値として
![{\displaystyle {\hat {\mu }}={\frac {1}{n}}\sum _{i=1}^{n}x_{i}}](https://wikimedia.org/api/rest_v1/media/math/render/svg/cee4ce53e038701485e1782671a780b6e2cdba8a)
![{\displaystyle {\hat {\sigma ^{2}}}={\frac {1}{n}}\sum _{i=1}^{n}(x_{i}-\mu )^{2}}](https://wikimedia.org/api/rest_v1/media/math/render/svg/cc1376f65e193b93120ffbd63847c29aa9424be2)
を得る。
ワイブル分布
Xi (i=1,..,n)が形状パラメータをβ、尺度パラメータをηとするワイブル分布に従うとする。このとき、対数尤度関数は
![{\displaystyle 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}^{n}x_{i}^{\beta }}](https://wikimedia.org/api/rest_v1/media/math/render/svg/94f0871f19225e40db411200440f565f921a4f5d)
であり、尤度方程式は
![{\displaystyle {\frac {\partial l(\eta ,\beta ,\mathbf {x} )}{\partial \eta }}=-{\frac {n\beta }{\eta }}-{\frac {\beta }{\eta ^{(\beta +1)}}}\sum _{i=1}^{n}x_{i}^{\beta }=0}](https://wikimedia.org/api/rest_v1/media/math/render/svg/71c10bde015e5d22c80976596a7e8544caf19384)
![{\displaystyle {\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}^{n}x_{i}^{\beta }+{\frac {1}{\eta ^{\beta }}}\sum _{i=1}^{n}\ln {x_{i}}x_{i}^{\beta }=0}](https://wikimedia.org/api/rest_v1/media/math/render/svg/476b44f800f30037ae37003849954a343fbdeb9b)
となる。 これらを整理すると最尤推定値ˆη、ˆβが満たすべき関係式
![{\displaystyle {\hat {\eta }}=\left({\frac {1}{n}}\sum _{i=1}^{n}x_{i}^{\hat {\beta }}\right)^{\frac {1}{\hat {\beta }}}}](https://wikimedia.org/api/rest_v1/media/math/render/svg/6079478f8864baad65942ce1eebe1aa8e218292c)
![{\displaystyle {\frac {1}{\hat {\beta }}}+{\frac {1}{n}}\sum _{i=1}^{n}\ln {x_{i}}-{\frac {\sum _{i=1}^{n}x_{i}^{\hat {\beta }}\ln {x_{i}}}{\sum _{i=1}^{n}x_{i}^{\hat {\beta }}}}=0}](https://wikimedia.org/api/rest_v1/media/math/render/svg/881f1a0f6c9c321dd1ebbff96f0d7c3ad6d731ce)
を得る。第二式を満たすˆβを数値的に求めれば、第一式よりˆηも定まる。
ガンマ分布
Xi (i=1,..,n)が形状パラメータをα、尺度パラメータをβとするガンマ分布に従うとする(X ∼ Γ(α, β))。このとき、対数尤度関数は
![{\displaystyle 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}^{n}x_{i}}](https://wikimedia.org/api/rest_v1/media/math/render/svg/4b83d05029c2ea92875e49dc9d3045ff712af497)
であり、尤度方程式は
![{\displaystyle {\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}](https://wikimedia.org/api/rest_v1/media/math/render/svg/bd312d25f7f9d673f6bc248b03d1a1dbc1910e11)
![{\displaystyle {\frac {\partial l(\alpha ,\beta ,\mathbf {x} )}{\partial \beta }}=-{\frac {n\alpha }{\beta }}+{\frac {1}{\beta ^{2}}}\sum _{i=1}^{n}x_{i}=0}](https://wikimedia.org/api/rest_v1/media/math/render/svg/bfba9b519ccff013570ed63341f35c4eb868eb7d)
となる。 ここではψ(α)はガンマ関数の対数微分であるディガンマ関数を表す。これらを整理すると最尤推定値ˆβ、ˆαが満たすべき関係式
![{\displaystyle {\hat {\beta }}={\frac {1}{\hat {\alpha }}}{\frac {1}{n}}\sum _{i=1}^{n}x_{i}}](https://wikimedia.org/api/rest_v1/media/math/render/svg/ec7c449dfda71ce64700f82464bdfe6d287a53e7)
![{\displaystyle {\hat {\alpha }}={\frac {{\frac {1}{n}}\sum _{i=1}^{n}x_{i}}{\left(\prod _{i=1}^{n}x_{i}\right)^{\frac {1}{n}}}}\exp {(\psi ({\hat {\alpha }}))}}](https://wikimedia.org/api/rest_v1/media/math/render/svg/9ee189f642296fd48f6691c693b6827db1c33341)
を得る。第二式を満たすˆαを数値的に求めれば、第一式よりˆβも定まる。
数値解法
尤度方程式が解析的に解けない場合、S(θ*)=0を満たすθ*∈ Θを数値的に求めることが必要となる。
ニュートン=ラフソン法
ニュートン=ラフソン法では、反復計算により、最適解θ*を求める。反復計算のkステップ目で求まったパラメータをθ(k)とする。スコア関数はテイラー展開により、
![{\displaystyle \mathbf {S} (\mathbf {x} ,{\boldsymbol {\theta }})\simeq \mathbf {S} (\mathbf {x} ,{\boldsymbol {\theta }}^{(k)})-I({\boldsymbol {\theta }}^{(k)})({\boldsymbol {\theta }}-{\boldsymbol {\theta }}^{(k)})}](https://wikimedia.org/api/rest_v1/media/math/render/svg/9d76d5e61b1b4c258776044050d686031474e04f)
と一次近似できる。ここでI(θ)は、
![{\displaystyle I({\boldsymbol {\theta }})=-{\frac {\partial ^{2}}{\partial {\boldsymbol {\theta }}\partial {\boldsymbol {\theta }}^{T}}}\ln {L({\boldsymbol {\theta }},\mathbf {x} )}}](https://wikimedia.org/api/rest_v1/media/math/render/svg/fbc79c2b9a174999a71958978b883e8656a3050f)
で与えられる、対数尤度関数のヘッセ行列の符号を変えた行列である。ニュートン=ラフソン法では、左辺をゼロとおくことで、θ(k+1)を与える更新式
![{\displaystyle {\boldsymbol {\theta }}^{(k+1)}={\boldsymbol {\theta }}^{(k)}+I({\boldsymbol {\theta }}^{(k)})^{-1}\mathbf {S} (\mathbf {x} ,{\boldsymbol {\theta }}^{(k)})}](https://wikimedia.org/api/rest_v1/media/math/render/svg/a7774cad04be35ee9d4da79d02c10a7bdda95536)
を定める。
ニュートン=ラフソン法は、最適解θ*の近傍で二次収束するため、収束が早い。すなわち、θ*の十分近くの適切な初期値を与えれば、
![{\displaystyle ||{\boldsymbol {\theta }}^{(k)}-{\boldsymbol {\theta }}^{\ast }||\leq K||{\boldsymbol {\theta }}^{(k)}-{\boldsymbol {\theta }}^{\ast }||^{2}}](https://wikimedia.org/api/rest_v1/media/math/render/svg/6edaaeb2bcf9cd45f2af89f59b5c26fd35f86ebe)
を満たす正の定数Kが存在する。
一方で、ニュートン=ラフソン法は各ステップで、対数尤度関数のヘッセ行列から定まるI(θ)の逆行列を計算する、もしくは、p次の連立方程式を解くことが必要となる。これらの計算量はO(p3)のオーダーであり、パラメータ数pが増えると、計算負荷が急激に増える。また、初期値の設定によっては、I(θ)は正定値とはならず、最適解θ*に収束しない場合がある。
フィッシャーのスコア法
ニュートン=ラフソン法においては、各ステップで負の対数尤度関数の二階微分であるI(θ)を計算する必要がある。このI(θ)を求める計算は、場合によっては煩雑となる。分布によっては、I(θ)の期待値であるフィッシャー情報行列
![{\displaystyle 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]}](https://wikimedia.org/api/rest_v1/media/math/render/svg/672a8f1946a064847286d4ace693b662f0d4c856)
が、より簡潔に求まるため、I(θ)をJ(θ)で代用し、反復計算を
![{\displaystyle {\boldsymbol {\theta }}^{(k+1)}={\boldsymbol {\theta }}^{(k)}+J({\boldsymbol {\theta }}^{(k)})^{-1}\mathbf {S} (\mathbf {x} ,{\boldsymbol {\theta }}^{(k)})}](https://wikimedia.org/api/rest_v1/media/math/render/svg/093cc1f53ff6241c135ae4a10b5d30eedd01e43c)
とする。この方法をフィッシャーのスコア法と呼ぶ。
フィッシャー情報行列は非負定値であるため、ニュートン=ラフソン法でのI(θ)の正定値性の問題を回避することができる。
脚注
参考文献
- Epps, T. W. (2013). Probability and Statistical Theory for Applied Researchers. World Scientific Pub Co Inc. ISBN 978-9814513159
- Lehmann, E. L. (1983). Theory of Point Estimation. John Wiley & Sons Inc. ISBN 978-0471058496
- Monahan, John F. (2011). Numerical Methods of Statistics. Cambridge Series in Statistical and Probabilistic Mathematics (2nd ed.). Cambridge University Press. ISBN 978-0521139519
関連項目