パラメータ推定(ニュートン・ラフソン法)
前回の続き。
「ニュートン・ラフソン法」を使って、2変数パラメータの最尤推定値を求めてみる。
以下のサイトを参考にさせていただいた。
http://www012.upp.so-net.ne.jp/doi/sas/numerical/Newton_Raphson/Newton_Raphson.pdf
正規分布の対数尤度関数を再掲。
データ
# 平均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とおいて、先ほどの式を計算した結果。
最尤値となるパラメータは、平均「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)
固有値分解の結果
最尤値をとるパラメータを与えたとき、ヘッセ行列の固有ベクトルは標準基底になっている。
また、各パラメータの2階偏微分の値が固有値になっている。
推定するパラメータが2つの場合、フィッシャー情報行列は以下として定義される。
今回は1セットのデータのみ用いており、パラメータ推定値も1セットのみであるから、フィッシャー情報行列は、対数尤度関数の2階偏微分にマイナスを掛けたものになる。
今回、各繰り返しにおけるフィッシャー情報量と対数尤度は以下のとおりであった。
平均を固定して、分散を変えてみると以下のとおりであった。
フィッシャー情報行列は、対数尤度関数のスコア関数(1階偏微分)の分散として定義されるから、以下のことがいえる。
- スコア関数の分散が大きい(=フィッシャー情報量が大きい)ほど、相対的に尤度は下がる
- スコア関数の分散が小さい(=フィッシャー情報量が小さい)ほど、相対的に尤度は上がる