Rで最尤法の復習 - データ分析メモと北欧生活

自分で尤度関数を書いて最尤法で推定する必要があったので、簡単に最尤法 (Maximum Likelihood Estimation)を復習。

読んで字のごとく、最も尤もらしいパラメーターを見つける推定法が最尤法である。

シンプルな生産関数の最尤推定

例えば、以下の生産関数のパラメータ(a,b,c)を推定したいとする。

$Y = aL^{b}K^{c}e$

Yは生産量、aは技術パラメータ、Lは投入労働量、Kは投入資本量、bとcはそれぞれパラメータで、eは誤差項である。
この誤差項は対数正規分布を取ると仮定する。

両辺の対数を取ると、

$\ln{Y} = \ln{a} + b\ln{a} + c\ln{K} + \varepsilon$ (1)

となる。 $e$ は対数正規分布を取るので、 $\varepsilon \equiv \ln{(e)}$ は正規分布 $N(0,\sigma^2)$ に従う。

正規分布の確率密度関数 $\phi(\varepsilon)$ は、

$\phi(\varepsilon) = \dfrac {1}{\sqrt{2\pi\sigma^2}}\exp\biggl\{ -\dfrac {(\varepsilon)^2}{2\sigma^2}\biggr\}$ (2)

である。

対数を取った生産関数(1)は $\varepsilon$ を左辺に置いて書くと、

$\varepsilon = \ln{Y} - \ln{a} - b\ln{a} - c\ln{K}$

なので、 $\phi(\varepsilon) = \phi( \ln{Y} - \ln{a} - b\ln{L} - c\ln{K})$ と書ける。

この確率密度関数は、パラメータ( $a,b,c, \sigma^{2}$ )を所与とした、Y, L, Kの関数だ。
つまり $a,b,c, \sigma^{2}$ がわかっていて、Y,L,Kが変化するという考え方。

$f(Y,L,K;a,b,c,\sigma^2) = \phi( \ln{Y} - \ln{a} - b\ln{L} - c\ln{K})$

しかし、データは既に取られていて手にある、つまりY,L,Kは固定。一方で、パラメータ(a,b,c, $\sigma^{2}$ )を探したい。
そこで、考え方としては確率密度関数を異なる見方、つまりデータを所与として、パラメータを求める関数という見方ができる。
これが尤度関数である。つまり固定されたデータがある状態で、パラメータを変化させて出た値はそのデータが現れる確率（尤もらしさ）である。このもっともらしさを最大化する、というのが最尤法の考え方である。

仮にデータは $i = 1,...,N$ の観測値があるとする。
誤差項 $e$ は、独立同分布(i.i.d)であると仮定して、同時確率密度関数は

$\prod_{i=1}^{N} f(Y_i,L_i,K_i;a,b,c,\sigma^2)$

対数を取ると、
$\ln(\prod_{i=1}^{N} f(Y_i,L_i,K_i;a,b,c,\sigma^2)) = \sum_{i=1}^{N} \ln{f(Y_i,L_i,K_i;a,b,c,\sigma^2)}$

これを対数尤度関数として定義する。

$\ln{Ln(\theta)} = \sum_{i=1}^{N} \ln{f(Y_i,L_i,K_i;\theta)} [tex: \theta$ はパラメータのベクトルを表している。簡便のため。 $\theta \equiv (a,b,c,\sigma^2)'$

この対数尤度関数を正規分布の確率密度関数(2)を使って特定化すると

$\ln{Ln(\theta)} = n\ln\frac{1}{\sqrt{2\pi\sigma^{2}}} - \sum_{i=1}^{n}\frac{(\ln{Y} - \ln{a} - b\ln{L} - c\ln{K})^{2}}{2\sigma^{2}}$

ここで、 $\ln\frac{1}{\sqrt{2\pi\sigma^{2}}} = \ln{(2\pi\sigma^2)^{-\frac{1}{2}}}$ なので、

$\ln{Ln(\theta)} = - \frac{n}{2}\ln{2\pi} - \frac{n}{2}\ln{2\sigma^{2}} - \sum_{i=1}^{n}\frac{(\ln{Y} - \ln{a} - b\ln{L} - c\ln{K})^{2}}{2\sigma^{2}}$

と書くことが出来る。第一項と第二項を分けたのは、第二項はパラメータ $\sigma^2$ を含めているが、第一項は定数項だからである。

データ(Y,L,K)が与えられている状態で、上記の対数尤度関数を最大にするようなパラメータが推定されるパラメータである。

Rコード

以下のRコードではdfというデータフレームが与えられていると仮定する。（最後にデータ生成過程を記載)

まず、対数尤度関数を、パラメータとデータを引数に取る関数として書く。
この時、一つ目の引数はベクトルにして、パラメータのベクトルを取る。ここではbeta[1]がa, beta[2]がb, beta[3]がc, beta[4]が $sigma^2$ である。残りの引数は各変数のデータベクトルを指定する。
ここで、対数尤度関数に-1をかけておく。

mloglik = function(beta,Y,L,K){
     n = length(Y)    # 観測値の数 (n) を取得
     sum( (log(Y)-beta[1]-beta[2]*log(L)-beta[3]*log(K))^2 )/(2*beta[4]^2) + 
     n/2*log(2*pi) + n*log(beta[4]) # 対数尤度関数
    }

対数尤度関数に-1をかけたのは、以下で使用する最適化が最小値を探す関数だからである。つまり、ある関数の負の値の最小化は、その関数の最大化である。
最適化のために使う関数はnlmである。オプションでhessianをTRUEにすることで、分散共分散行列を計算できる。(分散共分散行列の計算は別の機会に復習しようと思う。)

  mlem = nlm(mloglik,c(1,.5,.5,1),Y=df$Y,L=df$L,K=df$K, hessian = TRUE) # 対数尤度関数、パラメータの初期値、使用するデータを記述する。

> mlem
$minimum
[1] 1051.162

$estimate
[1] 0.7866748 0.2965670 0.7651674 1.9805737

$gradient
[1] -0.0002764864 -0.0009104042 -0.0004586127  0.0004539268

$hessian
              [,1]          [,2]         [,3]          [,4]
[1,] 127.464181787  395.78665110  499.7236374  -0.006084502
[2,] 395.786651097 1770.28630333 1509.1484784  -0.087547950
[3,] 499.723637404 1509.14847836 2430.4368253  -0.121001231
[4,]  -0.006084502   -0.08754795   -0.1210012 254.800541914

$code
[1] 1

$iterations
[1] 19

以上がRでの最尤法のステップである。

データ生成過程（付録)

以下のようにデータを生成した。
パラメータの真の値は $a = 3, b = 0.3, c = 0.7, \sigma^2 = 2$ である。

set.seed(3)
  
  N = 500;a = 3;b = 0.3;c = 0.7
  
  L = rlnorm(n = N, 3, 2)
  K = rlnorm(n = N, 4, 2)
  
  Y = exp(log(a) + b*log(L) + c*log(K) + rnorm(N,0,2))
  
  df = data.frame(Y,L,K)

元ネタは前回と同じくEconometrics in Rの14.3節。でも少し間違っていると思う(書いてある通りadditive errorだとRコードの尤度関数にならない）ので、この記事ではmultiplicative errorにした。