パラメータ推定(ニュートン・ラフソン法)

前回の続き。

ニュートン・ラフソン法」を使って、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)

固有値分解の結果
f:id:ichou1:20181029193007p:plain

最尤値をとるパラメータを与えたとき、ヘッセ行列の固有ベクトルは標準基底になっている。
また、各パラメータの2階偏微分の値が固有値になっている。


推定するパラメータが2つの場合、フィッシャー情報行列は以下として定義される。



\begin{eqnarray}
I(\mu, \sigma^2)
&=& - E(\nabla^2 (\mu, \sigma^2 ))
\end{eqnarray}

今回は1セットのデータのみ用いており、パラメータ推定値も1セットのみであるから、フィッシャー情報行列は、対数尤度関数の2階偏微分にマイナスを掛けたものになる。



\begin{eqnarray}
I(\mu, \sigma^2)
&=& - \nabla^2 (\mu, \sigma^2 )
\end{eqnarray}



今回、各繰り返しにおけるフィッシャー情報量と対数尤度は以下のとおりであった。
f:id:ichou1:20181029205949p:plain

平均を固定して、分散を変えてみると以下のとおりであった。
f:id:ichou1:20181029212531p:plain

フィッシャー情報行列は、対数尤度関数のスコア関数(1階偏微分)の分散として定義されるから、以下のことがいえる。

  • スコア関数の分散が大きい(=フィッシャー情報量が大きい)ほど、相対的に尤度は下がる
  • スコア関数の分散が小さい(=フィッシャー情報量が小さい)ほど、相対的に尤度は上がる