尤度方程式

**尤度方程式**（ゆうどほうていしき、）とは、統計学において、対数尤度関数の極値条件を与える方程式の事。統計的推定法の一つである最尤法において、尤度関数を最大化する最尤推定値を求める際に用いられる。

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

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

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

:
L(\boldsymbol{\theta} | \boldsymbol{d})
= \prod_{i=1}^n f(X=d_i|\boldsymbol{\theta})

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

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

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

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


この方程式を**尤度方程式**という。左辺の勾配ベクトル：

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


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

## 例
### 正規分布
が平均を、分散をとする正規分布に従うとする（）。このとき、対数尤度関数は
:
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


を得る。

### ワイブル分布
が形状パラメータを、尺度パラメータをとするワイブル分布に従うとする。このとき、対数尤度関数は

:
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


を得る。第二式を満たすを数値的に求めれば、第一式よりも定まる。

### ガンマ分布
が形状パラメータを、尺度パラメータをとするガンマ分布に従うとする（）。このとき、対数尤度関数は

:
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 \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}))}


を得る。第二式を満たすを数値的に求めれば、第一式よりも定まる。

## 数値解法
尤度方程式が解析的に解けない場合、を満たすを数値的に求めることが必要となる。

### ニュートン＝ラフソン法
ニュートン＝ラフソン法では、反復計算により、最適解を求める。反復計算のkステップ目で求まったパラメータをとする。スコア関数はテイラー展開により、

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


と一次近似できる。ここでは、

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


で与えられる、対数尤度関数のヘッセ行列の符号を変えた行列である。ニュートン＝ラフソン法では、左辺をゼロとおくことで、を与える更新式

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


を定める。

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

:
||\boldsymbol{\theta}^-\boldsymbol{\theta}^{\ast}||
\leq K ||\boldsymbol{\theta}^-\boldsymbol{\theta}^{\ast}||^2


を満たす正の定数が存在する。

一方で、ニュートン＝ラフソン法は各ステップで、対数尤度関数のヘッセ行列から定まるの逆行列を計算する、もしくは、次の連立方程式を解くことが必要となる。これらの計算量はのオーダーであり、パラメータ数が増えると、計算負荷が急激に増える。また、初期値の設定によっては、は正定値とはならず、最適解に収束しない場合がある。

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

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



が、より簡潔に求まるため、をで代用し、反復計算を

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


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

フィッシャー情報行列は非負定値であるため、ニュートン＝ラフソン法でのの正定値性の問題を回避することができる。

## 脚注

## 参考文献
* 
* 
* 

## 関連項目
* 最尤法