尤度方程式のソースを表示

'''尤度方程式'''（ゆうどほうていしき、{{lang-en-short|likelihood equation}}）とは、[[統計学]]において、対数[[尤度関数]]の[[極値]]条件を与える[[方程式]]の事{{sfn|Lehmann|1983| loc=&sect;6}}{{sfn|Epps|2013| loc=&sect;7}}。[[推計統計学|統計的推定法]]の一つである[[最尤法]]において、[[尤度関数]]を最大化する[[最尤推定値]]を求める際に用いられる。

== 概要 ==
[[独立同分布]]を満たす <math>n</math> 個の[[確率変数]] <math>\boldsymbol{D} = \{D_i \mid i \in \{1, .., n\} \}</math> とその観測値 <math>\boldsymbol{d} = \{d_i \mid i \in \{1, .., n\} \}</math> を定義する。すなわち真の分布から <math>n</math> 個の観測値（データ）が[[無作為抽出]]された状況を考える。

ここで[[確率密度関数]] <math>f(X|\boldsymbol{\theta})</math> に従う確率モデルを導入する。ここで <math>\boldsymbol{\theta} = {(\theta_1, .., \theta_p)}</math> は分布パラメータ群であり、パラメータ空間{{math|&Theta; &sub; '''R'''<sup>''p''</sup>}}に値を持つ。この確率モデルが <math>\boldsymbol{d}</math> を最も良く説明する <math>\boldsymbol{\theta}</math> を求めたい。ゆえに[[最尤推定]]をおこなう。

このとき[[独立同分布]]条件により、[[尤度関数]] <math>L(\boldsymbol{\theta} | \boldsymbol{d})</math> と対数尤度関数 <math>l(\boldsymbol{\theta} | \boldsymbol{d})</math> は以下で定義される。

:<math>
L(\boldsymbol{\theta} | \boldsymbol{d})
= \prod_{i=1}^n f(X=d_i|\boldsymbol{\theta})
</math>
:<math>
l(\boldsymbol{\theta} | \boldsymbol{d})
= \ln{L(\boldsymbol{\theta}| \boldsymbol{d})}
=\sum_{i=1}^n \ln{f(X=d_i|\boldsymbol{\theta})}</math>

すなわちあるデータ群に対するモデルの尤度関数は、各観測値に対する尤度関数の積（対数尤度の場合は和）となる。

最尤法では対数尤度関数を最大化する <math>\boldsymbol{\theta}</math> が最尤推定値 <math>\hat{{\boldsymbol{\theta}}}</math> として定まる。このとき <math>\hat{{\boldsymbol{\theta}}}</math> は次の極値条件を満たす。

:<math>
\frac{\partial}{\partial \boldsymbol{\theta}}
l(\boldsymbol{\theta} | \boldsymbol{d})
=\mathbf{0}
</math>

この方程式を'''尤度方程式'''という。左辺の[[勾配 (ベクトル解析)|勾配ベクトル]]：

:<math>
\mathbf{S}( \boldsymbol{d},\boldsymbol{\theta})
:=\frac{\partial}{\partial \boldsymbol{\theta}}
l(\boldsymbol{\theta} | \boldsymbol{d})
</math>

は、'''スコア関数'''、もしくは単に'''スコア'''と呼ばれる。多くの場合、最尤推定値の推定は、尤度方程式を解く問題、すなわち、スコアをゼロとするパラメータ{{math|'''&theta;'''&isin; &Theta;}}を求める問題に帰着する。

== 例 ==
=== 正規分布 ===
{{math|''X<sub>i</sub>'' (''i''{{=}}1,..,''n'')}}が平均を{{mvar|&mu;}}、分散を{{mvar|&sigma;<sup>2</sup>}}とする[[正規分布]]に従うとする（{{math|''X'' &sim; N(&mu;, &sigma;<sup>2</sup>)}}）。このとき、対数尤度関数は
:<math>
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
</math>

であり、尤度方程式は

:<math>
\frac{\partial l(\mu,\sigma^2,\mathbf{x})}{\partial \mu}=\frac{1}{\sigma^2}\sum_{i=1}^n(x_i-\mu)=0
</math>
:<math>
\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
</math>

となる。これらを整理すると最尤推定値として

:<math>
\hat{\mu}=\frac{1}{n} \sum_{i=1}^nx_i
</math>
:<math>
\hat{\sigma^2}=\frac{1}{n}\sum_{i=1}^n(x_i-\mu)^2
</math>

を得る。

=== ワイブル分布 ===
{{math|''X<sub>i</sub>'' (''i''{{=}}1,..,''n'')}}が形状パラメータを{{mvar|&beta;}}、尺度パラメータを{{mvar|&eta;}}とする[[ワイブル分布]]に従うとする。このとき、対数尤度関数は

:<math>
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}
</math>

であり、尤度方程式は

:<math>
\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
</math>
:<math>
\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
</math>

となる。
これらを整理すると最尤推定値{{math|{{hat|&eta;}}}}、{{math|{{hat|&beta;}}}}が満たすべき関係式

:<math>
\hat{\eta}=\left( \frac{1}{n}\sum_{i=1}^n x_i^{\hat{\beta}} \right)^{\frac{1}{\hat{\beta}}}
</math>
:<math>
\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
</math>

を得る。第二式を満たす{{math|{{hat|&beta;}}}}を数値的に求めれば、第一式より{{math|{{hat|&eta;}}}}も定まる。

=== ガンマ分布 ===
{{math|''X<sub>i</sub>'' (''i''{{=}}1,..,''n'')}}が形状パラメータを{{mvar|&alpha;}}、尺度パラメータを{{mvar|&beta;}}とする[[ガンマ分布]]に従うとする（{{math|''X'' &sim; &Gamma;(&alpha;, &beta;)}}）。このとき、対数尤度関数は

:<math>
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
</math>

であり、尤度方程式は

:<math>
\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
</math>
:<math>
\frac{\partial l(\alpha,\beta,\mathbf{x})}{\partial \beta}=
-\frac{n \alpha}{\beta}+\frac{1}{\beta^2}\sum_{i=1}^nx_i=0
</math>

となる。
ここでは{{math|&psi;(&alpha;)}}はガンマ関数の[[対数微分]]である[[ディガンマ関数]]を表す。これらを整理すると最尤推定値{{math|{{hat|&beta;}}}}、{{math|{{hat|&alpha;}}}}が満たすべき関係式

:<math>
\hat{\beta}=\frac{1}\hat{{\alpha}} \frac{1}{n}\sum_{i=1}^nx_i
</math>
:<math>
\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}))}
</math>

を得る。第二式を満たす{{math|{{hat|&alpha;}}}}を数値的に求めれば、第一式より{{math|{{hat|&beta;}}}}も定まる。

== 数値解法 ==
尤度方程式が解析的に解けない場合、{{math|'''S'''('''&theta;*'''){{=}}'''0'''}}を満たす{{math|'''&theta;*'''&isin; &Theta;}}を数値的に求めることが必要となる{{sfn|Monahan|2011| loc=&sect;9}}。

=== ニュートン＝ラフソン法 ===
[[ニュートン＝ラフソン法]]では、反復計算により、最適解{{math|'''&theta;*'''}}を求める。反復計算のkステップ目で求まったパラメータを{{math|'''&theta;'''<sup>(''k'')</sup>}}とする。スコア関数は[[テイラー展開]]により、

:<math>
\mathbf{S}(\mathbf{x},\boldsymbol{\theta})\simeq 
\mathbf{S}(\mathbf{x},\boldsymbol{\theta}^{{(k)}})-I(\boldsymbol{\theta}^{(k)})(\boldsymbol{\theta}-\boldsymbol{\theta}^{(k)})
</math>

と[[一次近似]]できる。ここで{{math|''I''('''&theta;''')}}は、

:<math>
I(\boldsymbol{\theta})=
-\frac{\partial^2}{\partial \boldsymbol{\theta}\partial \boldsymbol{\theta}^T}\ln{L(\boldsymbol{\theta}, \mathbf{x})}
</math>

で与えられる、対数尤度関数の[[ヘッセ行列]]の符号を変えた行列である。ニュートン＝ラフソン法では、左辺をゼロとおくことで、{{math|'''&theta;'''<sup>(''k''+1)</sup>}}を与える更新式

:<math>
\boldsymbol{\theta}^{(k+1)}
=\boldsymbol{\theta}^{(k)}+
I(\boldsymbol{\theta}^{(k)})^{-1}\mathbf{S}(\mathbf{x},\boldsymbol{\theta}^{{(k)}})
</math>

を定める。

ニュートン＝ラフソン法は、最適解{{math|'''&theta;*'''}}の近傍で二次収束するため、収束が早い。すなわち、{{math|'''&theta;*'''}}の十分近くの適切な初期値を与えれば、

:<math>
||\boldsymbol{\theta}^{{(k)}}-\boldsymbol{\theta}^{\ast}||
\leq K ||\boldsymbol{\theta}^{{(k)}}-\boldsymbol{\theta}^{\ast}||^2
</math>

を満たす正の定数{{mvar|''K''}}が存在する。

一方で、ニュートン＝ラフソン法は各ステップで、対数尤度関数のヘッセ行列から定まる{{math|''I''('''&theta;''')}}の逆行列を計算する、もしくは、{{mvar|p}}次の連立方程式を解くことが必要となる。これらの[[計算量]]は{{math|''O''(''p''<sup>3</sup>)}}の[[ランダウの記号|オーダー]]であり、パラメータ数{{mvar|p}}が増えると、計算負荷が急激に増える。また、初期値の設定によっては、{{math|''I''('''&theta;''')}}は[[行列の定値性|正定値]]とはならず、最適解{{math|'''&theta;*'''}}に収束しない場合がある。

=== フィッシャーのスコア法 ===
ニュートン＝ラフソン法においては、各ステップで負の対数尤度関数の二階微分である{{math|''I''('''&theta;''')}}を計算する必要がある。この{{math|''I''('''&theta;''')}}を求める計算は、場合によっては煩雑となる。分布によっては、{{math|''I''('''&theta;''')}}の[[期待値]]である[[フィッシャー情報行列]]

:<math>
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 ]

</math>

が、より簡潔に求まるため、{{math|''I''('''&theta;''')}}を{{math|''J''('''&theta;''')}}で代用し、反復計算を

:<math>
\boldsymbol{\theta}^{(k+1)}
=\boldsymbol{\theta}^{(k)}+
J(\boldsymbol{\theta}^{(k)})^{-1}\mathbf{S}(\mathbf{x},\boldsymbol{\theta}^{{(k)}})
</math>

とする。この方法を'''フィッシャーのスコア法'''と呼ぶ。

フィッシャー情報行列は非負定値であるため、ニュートン＝ラフソン法での{{math|''I''('''&theta;''')}}の正定値性の問題を回避することができる。

== 脚注 ==
{{Reflist}}

== 参考文献 ==
* {{cite book
|last1      = Epps 
|first1     = T. W.
|year       = 2013
|title      = Probability and Statistical Theory for Applied Researchers 
|publisher  = World Scientific Pub Co Inc
|isbn       = 978-9814513159
|ref        = harv
}}
* {{cite book
|last1      = Lehmann
|first1     = E. L.
|year       = 1983
|title      = Theory of Point Estimation
|publisher  = John Wiley & Sons Inc
|isbn       =  978-0471058496
|ref        = harv
}}
* {{cite book
|last1      = Monahan 
|first1     = John F.
|year       = 2011
|title      = Numerical Methods of Statistics
|series     = Cambridge Series in Statistical and Probabilistic Mathematics
|edition    = 2nd
|publisher  = Cambridge University Press
|isbn       = 978-0521139519
|ref        = harv
}}

== 関連項目 ==
* [[最尤法]]

{{統計学}}
{{DEFAULTSORT:ゆうとほうていしき}}
[[Category:推計統計学]]
[[Category:方程式]]
[[Category:数学に関する記事]]