パラメータ推定（ニュートン・ラフソン法）

前回の続き。

「ニュートン・ラフソン法」を使って、2変数パラメータの最尤推定値を求めてみる。
以下のサイトを参考にさせていただいた。
http://www012.upp.so-net.ne.jp/doi/sas/numerical/Newton_Raphson/Newton_Raphson.pdf

正規分布の対数尤度関数を再掲。

対数尤度関数（正規分布）

$\begin{eqnarray} l(\mu, \sigma^2 ) &=& \log L(\mu, \sigma^2 ) \\ &=& \log \displaystyle \prod_{ i = 1 }^n f(x_i | \mu, \sigma^2 ) \\ &= & - \frac{n}{2} \log 2\pi - \frac{n}{2} \log \sigma^2 - \frac{1}{2\sigma^2} \displaystyle \sum_{ i = 1 }^{ n } (x_i - \mu)^2 \end{eqnarray}$

対数尤度関数をパラメータで1階偏微分する。

$\begin{eqnarray} \nabla l(\mu, \sigma^2) &=& \left( \begin{array}{c} \frac{ \partial l(\mu, \sigma^2 ) }{ \partial \mu } \\ \frac{ \partial l(\mu, \sigma^2 ) }{ \partial \sigma^2 } \end{array} \right) \\ &=& \left( \begin{array}{c} \frac{1}{\sigma^2} \displaystyle \sum_{ i = 1 }^{ n } (x_i - \mu) \\ - \frac{n}{2\sigma^2} + \frac{1}{2(\sigma^2)^2} \displaystyle \sum_{ i = 1 }^{ n } (x_i - \mu)^2 \end{array} \right) \end{eqnarray}$

2階偏微分を求める。

$\begin{eqnarray} \nabla^2 l(\mu, \sigma^2) &=& \left( \begin{array}{cc} \frac{ \partial }{ \partial \mu } \frac{ \partial l(\mu, \sigma^2 )}{ \partial \mu } && \frac{ \partial }{ \partial \sigma^2 } \frac{ \partial l(\mu, \sigma^2 )}{ \partial \mu } \\ \frac{ \partial }{ \partial \mu } \frac{ \partial l(\mu, \sigma^2 )}{ \partial \sigma^2 } && \frac{ \partial }{ \partial \sigma^2 } \frac{ \partial l(\mu, \sigma^2 )}{ \partial \sigma^2 } \end{array} \right) \\ &=& \left( \begin{array}{cc} - \frac{n}{\sigma^2} && - \frac{1}{(\sigma^2)^2} \displaystyle \sum_{ i = 1 }^{ n } (x_i - \mu) \\ - \frac{1}{(\sigma^2)^2} \displaystyle \sum_{ i = 1 }^{ n } (x_i - \mu) && \frac{n}{2(\sigma^2)^2 } - \frac {1} {(\sigma^2)^3} \displaystyle \sum_{ i = 1 }^{ n } (x_i - \mu)^2 \end{array} \right) \end{eqnarray}$

繰り返し計算の回数をkとおくと、漸化式は以下のとおり。

$\begin{eqnarray} \left( \begin{array}{c} \mu_{k+1} \\ \sigma_{k+1}^2 \end{array} \right) &=& \left( \begin{array}{c} \mu_k \\ \sigma^2_k \end{array} \right) - \left( \nabla^2 l(\mu_k, \sigma_k^2 ) \right)^{-1} \nabla l(\mu_k, \sigma_k^2 ) \\ &=& \left( \begin{array}{c} \mu_k \\ \sigma^2_k \end{array} \right) - \frac{\nabla l(\mu_k, \sigma_k^2 )}{\nabla^2 l(\mu_k, \sigma_k^2 ) } \end{eqnarray}$

検算。

データ

# 平均2、分散5の正規分布に従うサンプルを200個
set.seed(123); data <- rnorm(200, mean=2, sd=sqrt(5))

  [1]  0.746738355  1.485307487  5.485377748  2.157661556  2.289096164
  [6]  5.835001897  3.030639969 -0.828762916  0.464150333  1.003469540
 [11]  4.737130109  2.804568177  2.896152207  2.247494077  0.757101438
(snip)
[186]  1.559101297  4.481857217  2.189478346  3.686115522  0.883549109
[191]  2.479514290  1.273980231  2.211495199 -0.002093333 -0.931041334
[196]  6.465904894  3.343225764 -0.797927823  0.633391465 -0.650814055

初期パラメータは、平均1、分散1とおいて、先ほどの式を計算した結果。

f:id:ichou1:20181029185002p:plain
最尤値となるパラメータは、平均「1.980836」、分散「4.425514」

検算。Rの関数で計算した結果。

# 準ニュートン法による最適解
> optim(c(1, 1), fn=log_likelihood_norm, obs=data, method="BFGS", control=list(fnscale=-1), hessian=TRUE)

結果

$par
[1] 1.980839 4.425508

$value
[1] -432.5263

$counts
function gradient 
      31       10 

$convergence
[1] 0

$message
NULL

$hessian
              [,1]          [,2]
[1,] -4.519255e+01  3.119993e-05
[2,]  3.119993e-05 -5.105931e+00

小数点第五位まで一致していることが確認できた。

また、k回目の最尤値でのヘッセ行列を固有値分解してみる。

H_0 <- matrix(c(-200.000000 , -196.167181 , -196.167181 , -977.510617  ), nrow=2)
H_1 <- matrix(c(-152.232794 ,  -35.661895 ,  -35.661895 , -340.743870  ), nrow=2)
H_2 <- matrix(c(-108.715893 ,   -7.280127 ,   -7.280127 , -113.101423  ), nrow=2)
H_3 <- matrix(c( -79.549506 ,   -1.429134 ,   -1.429134 ,  -39.900363  ), nrow=2)
H_4 <- matrix(c( -61.130669 ,   -0.254284 ,   -0.254284 ,  -15.933080  ), nrow=2)
H_5 <- matrix(c( -50.656386 ,   -0.036104 ,   -0.036104 ,   -7.966423  ), nrow=2)
H_6 <- matrix(c( -46.162098 ,   -0.002919 ,   -0.002919 ,   -5.555946  ), nrow=2)
H_7 <- matrix(c( -45.231581 ,   -0.000058 ,   -0.000058 ,   -5.123588  ), nrow=2)
H_8 <- matrix(c( -45.192558 ,   -0.000000 ,   -0.000000 ,   -5.105933  ), nrow=2)
H_9 <- matrix(c( -45.192490 ,   -0.000000 ,   -0.000000 ,   -5.105903  ), nrow=2)