計量分析で最初に習うOLSは性質がよく計算しやすい推定ですが、非線形モデルの推定に使えません。非線形関数を推定する場合は、その非線形関数を推定するパッケージを用いない限りは、最尤法が手っ取り早い選択になります¹。最尤法は、一般化線形回帰モデル（GLM）やベイズ統計学の理解の前提になるなど、統計学の学習にも避けて通ることができない手法です。

実用上は尤度関数を書いて、何らかの凸最適化のアルゴリズムを実装して走らせるか、nlmかoptimといったパラメーターの最適化をしてくれる関数を呼べば簡単に最適点と最適点におけるヘッシアンを推定でき、ヘッシアンの意味を理解していたら標準誤差を計算することができ、それで推定したパラメーターを検定することができます。理解のために実際に、最尤法を試してみましょう。

1 対数尤度関数

最尤法は、観測値それぞれが抽出される確率（密度）を乗じた同時確率（密度）を尤度と定義し、尤度が最大化になるように確率モデルのパラメーターを推定する方法です。ただし、パラメーターの標準誤差を計算する都合などで、尤度ではなく対数尤度関数を用いて推定を行ないます。

例えば平均\(\mu\)，分散\(\sigma^2\)の正規分布を推定するときの尤度関数\(\mathfrak{L}\)は、観測値を\(y\)，観測数を\(n\)として、

\[ \begin{aligned} \mathfrak{L}(\mu, \sigma^2) &= \frac{1}{\sqrt{2 \pi \sigma^2}}\exp \big( - \frac{(y_1 - \mu)^2}{2\sigma^2} \big ) \cdots \frac{1}{\sqrt{2 \pi \sigma^2}}\exp \big( -\frac{(y_n - \mu)^2}{2\sigma^2} \big) \\ &= \big( \frac{1}{\sqrt{2\pi}} \big)^n \big( \frac{1}{\sqrt{\sigma^2}} \big)^n \exp{ \big( - \frac{1}{2\sigma^2} \sum^n_{i=1} (y_i-\mu)^2 \big) } \end{aligned} \]

となり、これを対数化したものを最適化の目的関数とします。

\[ \begin{aligned} &\log{ \mathfrak{L}(\mu, \sigma^2)} \\ &= \log{ (2\pi)^{-\frac{n}{2}}} + \log{ (\sigma^2)^{-\frac{n}{2}}} + \log{ \exp \big( -\frac{1}{2\sigma^2} \sum^n_{i=1} (y_i - \mu)^2 \big) } \\ &= - \frac{n}{2} \log{2\pi} - \frac{n}{2} \log{\sigma^2} - \frac{1}{2\sigma^2} \sum^n_{i=1} (y_i - \mu)^2 \end{aligned} \]

2 対数尤度関数のグラディエント

平均\(\mu\)，分散\(\sigma^2\)の最大化の一階条件を整理しましょう。

\[ \begin{aligned} \frac{ \partial \log{ \mathfrak{L}(\mu, \sigma^2)} }{ \partial \mu } &= - (-1) \cdot 2 \cdot \frac{1}{2\sigma^2} \sum^n_{i=1} (y_i - \mu) = \frac{1}{\sigma^2} \sum^n_{i=1} (y_i - \mu) = 0 \\ \frac{ \partial \log{ \mathfrak{L}(\mu, \sigma^2)} }{ \partial \sigma^2 } &= -\frac{n}{2\sigma^2} + \frac{1}{2 ( \sigma^2 )^2 } \sum^n_{i=1} (y_i - \mu)^2 = 0 \end{aligned} \]

この連立方程式の解が求めるものです。

3 ニュートン・ラフソン法で推定

実践的ではないですが、理解のために古典的なニュートン・ラフソン法で連立方程式を解いてみましょう。このアルゴリズムは手順が簡単で収束が速いと言う利点があります。

3.1 解析学的な意味

\(\mu\)と\(\sigma^2\)の2元の連立方程式をベクトル値関数\(g(\mu, \sigma^2)\)にまとめます。

\[ \begin{aligned} g(\mu, \sigma^2) &= \begin{bmatrix} \frac{1}{\sigma^2} \sum^n_{i=1} (y_i - \mu) \\ -\frac{n}{2\sigma^2} + \frac{1}{2 ( \sigma^2 )^2 } \sum^n_{i=1} (y_i - \mu)^2 \end{bmatrix} = 0 \end{aligned} \]

テイラー展開で\(f(\mu, \sigma^2)\)の一次近似式をつくると

\[ \begin{aligned} g(\mu, \sigma^2) = g(\mu_0, \sigma_0^2) + \frac{\partial g}{\partial(\mu_0, \sigma_0^2)} \begin{bmatrix} \mu_0 - \mu \\ \sigma_0^2 - \sigma^2 \end{bmatrix} \end{aligned} \]

になり、この零点を整理すると以下になります。

\[ \begin{aligned} \begin{bmatrix} \mu \\ \sigma^2 \end{bmatrix} = \begin{bmatrix} \mu_0 \\ \sigma_0^2 \end{bmatrix} - \frac{\partial g}{\partial(\mu_0, \sigma_0^2)}^{-1} g(\mu_0, \sigma_0^2) \end{aligned} \]

\(g\)が凸関数であれば、右辺の値は\(g(\mu_0, \sigma_0^2)\)よりも\(g(\mu, \sigma^2)\)に近づくので、右辺の計算を繰り返し行なえば解になる不動点に辿り着けます。

3.2 実際に計算するコード

上述の説明をだいたい忠実に実装したコードは以下になります。

# サンプル
set.seed(31)
y <- rnorm(30, mean=12, sd=7)
n <- length(y)
# 連立方程式を設定
g1 <- expression(sum(y-mu)/s2)
g2 <- expression(-n/(2*s2) + sum((y-mu)^2)/(2*(s2^2)))
# muとs2で一階微分を作っておく
h11 <- expression(-n/s2)
h12 <- expression(-sum(y - mu)/(s2^2))
h21 <- expression(-sum(y - mu)/(s2^2))
h22 <- expression(n/(2*(s2^2)) - sum( (y - mu)^2 )/(s2^3))

# 初期値を設定
mu <- 11 # 正規分布の平均値
s2 <- 16 # 正規分布の分散
for(i in 1:10){
    # (mu, s2)を2x1行列にする
    m <- matrix(c(mu, s2), 2, 1)
    # (mu, s2)で評価したg1とg2の値を2x1行列にする
    g <- matrix(c(eval(g1), eval(g2)), 2, 1)
    # ヤコビアンに(mu, s2)を代入した行列を作る
    j <- matrix(c(eval(h11), eval(h21), eval(h12), eval(h22)), 2, 2)
    # 行列mから、ヤコビアンの逆行列と評価式を乗じたものを引く
    m <- m - solve(j) %*% g
    # 行列mを(x, y)に展開しておく
    mu <- m[1]; s2 <- m[2]
    print(sprintf("[%d] (mu,s2)=(%f,%f)", i, mu, s2))
}

[1] "[1] (mu,s2)=(11.666249,21.695600)"
[1] "[2] (mu,s2)=(11.920991,28.383352)"
[1] "[3] (mu,s2)=(12.009620,34.606119)"
[1] "[4] (mu,s2)=(12.031770,38.413427)"
[1] "[5] (mu,s2)=(12.034435,39.431333)"
[1] "[6] (mu,s2)=(12.034508,39.488142)"
[1] "[7] (mu,s2)=(12.034508,39.488307)"
[1] "[8] (mu,s2)=(12.034508,39.488307)"
[1] "[9] (mu,s2)=(12.034508,39.488307)"
[1] "[10] (mu,s2)=(12.034508,39.488307)"

7順目ぐらいで収束して解が得られています。最適化までの経路をプロットすると以下になります。横軸が\(\sigma^2\)ではなくて\(\sigma\)なのに注意してください。

4 パラメーターの標準誤差

実際のところ分析では色々と検定する事になるので、パラメーターの標準誤差の算出方法を確認しましょう。

4.1 導出

テキストに沿い²、対数尤度関数のモーメントの性質を用いて、パラメーターの推定量が漸近的に正規分布に従うことと、その分散共分散行列を導出します。

4.1.1 対数尤度関数のモーメントの性質

まず、2つのパラメーター\(\mu, \sigma\)をベクトル\(\theta\)に置きなおします。

\[ \theta = \begin{bmatrix} \mu \\ \sigma^2 \end{bmatrix} \]

対数尤度関数を微分した式であるグラディエントを以下のように置きます。グラディエントはベクトルです。

\[ g_i = \frac{\partial}{\partial \theta} \ln f(y_i|\theta) \]

グラディエントの式をそれぞれ微分したヘッシアンを以下のように置きます。ヘッシアンは行列です。

\[ H_i = \frac{\partial^2}{\partial \theta \partial {}^t\!\theta} \ln f(y_i|\theta) \]

観察できない真のパラメーターを\(\theta_0\)と置きます。

次に、ルベーグの収束定理から明らかな気もちょっとするのですが、積分と微分の順序が交換できることを確認します。観測値の下限を\(A(\theta_0)\)、上限を\(B(\theta_0)\)と便宜的に置きます。尤度関数は確率密度関数なので、

\[ \int_{A(\theta_0)}^{B(\theta_0)} f(y_i|\theta) dy_i = 1 \]

となります。これをLeibniz integral ruleを用いて微分します。

\[\begin{align} &\frac{\partial}{\partial \theta_0} \int_{A(\theta_0)}^{B(\theta_0)} f(y_i|\theta) dy_i = 0\\ &= \int_{A(\theta_0)}^{B(\theta_0)} \frac{\partial}{\partial \theta_0} f(y_i|\theta) dy_i \\ &+ f(B(\theta_0)|\theta_0) \frac{\partial B(\theta_0)}{\partial \theta_0} - f(A(\theta_0)|\theta_0) \frac{\partial A(\theta_0)}{\partial \theta_0} \end{align}\]

パラメーターに観測値の上限と下限が影響されない場合は、第2項と第3項はゼロになります。一様分布は非ゼロですが、よく用いられる尤度関数はゼロになります。正則条件と呼ばれます。これで、微分と積分の順番が入れ替え可能だとわかります。

条件を満たすとき積分の上限と下限は気にしなくてよくなるので、\(A(\theta_0)\)と\(B(\theta_0)\)は省略します。

\[\begin{align} &\int \frac{\partial}{\partial \theta_0} f(y_i|\theta_0) dy_i \\ &= \int \bigg( \frac{\partial}{\partial \theta_0} \ln f(y_i|\theta_0) \bigg) f(y_i|\theta_0) dy_i \\ &= E\bigg[ \frac{\partial}{\partial \theta_0} \ln f(y_i|\theta_0) \bigg] = 0 \end{align}\]

対数微分の性質と、期待値の定義を使いました。

二階微分も微分と積分の順番が入れ替え可能なので、整理整頓ができます。

\[\begin{align} &\frac{\partial^2}{\partial \theta_0 \partial {}^t\!\theta_0} \int f(y_i|\theta_0) dy_i\\ &= \int \bigg[ \frac{\partial^2}{\partial \theta_0 \partial {}^t\!\theta_0} f(y_i|\theta_0) \bigg] dy_i \\ &= \int \bigg( \frac{\partial}{\partial {}^t\!\theta_0} \bigg[ \bigg( \frac{\partial}{\partial \theta_0} \ln f(y_i|\theta_0) \bigg) f(y_i|\theta_0) \bigg] dy_i \\ &= \int \bigg[ \bigg( \frac{\partial^2}{\partial \theta_0 \partial {}^t\! \theta_0} \ln f(y_i|\theta_0) \bigg) f(y_i|\theta_0) + \frac{\partial}{\partial \theta_0} \ln f(y_i|\theta_0) \frac{\partial}{\partial \theta_0} f(y_i|\theta_0) \bigg] dy_i \\ &= \int \bigg[ \bigg( \frac{\partial^2}{\partial \theta_0 \partial {}^t\! \theta_0} \ln f(y_i|\theta_0) \bigg) f(y_i|\theta_0) + \frac{\partial}{\partial \theta_0} \ln f(y_i|\theta_0) \frac{\partial}{\partial \theta_0} f(y_i|\theta_0) \bigg] dy_i \\ &= 0 \end{align}\]

積分の中の第一項の積分と、第二項の積分は等しいので、統合で結びます。

\[\begin{align} -\int \bigg( \frac{\partial^2}{\partial \theta_0 \partial {}^t\! \theta_0} \ln f(y_i|\theta_0) \bigg) f(y_i|\theta_0) dy_i &= \int \frac{\partial}{\partial \theta_0} \ln f(y_i|\theta_0) \frac{\partial}{\partial \theta_0} f(y_i|\theta_0) dy_i \\ -\int \bigg( \frac{\partial^2}{\partial \theta_0 \partial {}^t\! \theta_0} \ln f(y_i|\theta_0) \bigg) f(y_i|\theta_0) dy_i &= \int \frac{\partial}{\partial \theta_0} \ln f(y_i|\theta_0) \bigg( \frac{\partial}{\partial \theta_0} \ln f(y_i|\theta_0) \bigg) f(y_i|\theta_0) dy_i\\ -E\bigg[ \frac{\partial^2}{\partial \theta_0 \partial {}^t\! \theta_0} \ln f(y_i|\theta_0) \bigg] &= E\bigg[\frac{\partial}{\partial \theta_0} \ln f(y_i|\theta_0) \frac{\partial}{\partial \theta_0} \ln f(y_i|\theta_0) \bigg] \end{align}\]

もう一度、対数微分の性質と、期待値の定義を使いました。左辺は対数尤度関数のヘッシアンの符号を変えたものの期待値、右辺はグラディエントの分散になり、この二つは等しくなります。

\[ -E[H(y_i|\theta_0)] = VAR[g(y_i|\theta_0)] \]

4.1.2 パラメーターの分散共分散行列

上のコード例でのヤコビアン\(\partial g/\partial \theta\)は対数尤度関数を2階偏微分した行列で、対数尤度関数のヘッシアン\(H(\theta)\)になるのに注意して、真のパラメーター\(\theta_p\)の周りで、テイラーの定理を用います。2階微分の引数\(\bar{\theta}\)は、平均値の定理から存在することが分かる、\(\hat{\theta}\)と\(\theta_p\)の間にある値です。

\[ g(\hat{\theta}) = g(\theta_p) + H(\bar{\theta})(\hat{\theta} - \theta_p) = 0 \]

\[ \hat{\theta} - \theta_p = - H(\bar{\theta})^{-1} g(\theta_p) \]

両辺に\(\sqrt{n}\)を乗じます。

\[ \sqrt{n} (\hat{\theta} - \theta_p) = - H(\bar{\theta})^{-1} \sqrt{n} g(\theta_p) \DeclareMathOperator{\plim}{plim} \]

\(\plim \hat{\theta} - \theta_p = 0\)と言うことは、\(\hat{\theta}\)と\(\theta_p\)の間にある\(\bar{\theta}\)も\(\theta_p\)に確率収束するので、\(\sqrt{n} (\hat{\theta} - \theta_p)\)が\(- H(\theta_p)^{-1} \sqrt{n} g(\theta_p)\)に分布収束することが分かります。

\[ \sqrt{n} (\hat{\theta} - \theta_p) \xrightarrow{d} - H(\theta_p)^{-1} \sqrt{n} g(\theta_p) \]

\(g(\theta_p)\)を\(n\)で割り、\(H(\theta_p)^{-1}\)に\(n\)を乗じます。\(\bar{g}(\theta_p)=g(\theta_p)/n\)に注意すると、

\[ \sqrt{n} (\hat{\theta} - \theta_p) \xrightarrow{d} - \big[ \frac{1}{n} H(\theta_p) \big]^{-1} \sqrt{n} \bar{g}(\theta_p) \]

と書けます。

\(\sqrt{n}\bar{g}(\theta_p)\)は乱数の平均値になるので、Lindberg–Lévy中心極限定理を用いることができ、

\[ \sqrt{n} \bar{g}(\theta_p) \xrightarrow{d} \mathcal{N} \bigg \{ 0, -E \big[ \frac{1}{n} H(\theta_p) \big] \bigg \} \]

と書けます。ヘッシアンの期待値の負の値がグラディエントの分散になる、対数尤度関数のモーメントの性質を用いました。

大数の弱法則から\(\plim (-1/n)H(\theta_p) = - E [ (-1/n)H(\theta_p) ]\)なので、

\[ \begin{aligned} &- \big[ \frac{1}{n} H(\theta_p) \big]^{-1} \sqrt{n} \bar{g}(\theta_p) \xrightarrow{d} \\ &\mathcal{N} \bigg \{ 0, \bigg [ -E \big[ \frac{1}{n} H(\theta_p) \big] \bigg ]^{-1} \bigg [ -E \big[ \frac{1}{n} H(\theta_p) \big] \bigg ] \bigg [ -E \big[ \frac{1}{n} H(\theta_p) \big ] \bigg ]^{-1} \bigg \} \end{aligned} \]

になり、

\[ \sqrt{n} (\hat{\theta} - \theta_p) \xrightarrow{d} \mathcal{N} \bigg \{ 0, -E \big[ \frac{1}{n} H(\theta_p) \big]^{-1} \bigg \} \]

パラメーターの真の値と推定量の差の分散が、ヘッシアンの逆行列の対角成分に漸近することが分かります。

4.2 実際に計算するコード

説明は長いわけですが、計算は簡単です。上のコードのヤコビアンjがヘッシアンになるので、分散共分散行列は以下のように計算できます。

vcov <- -solve(j)
vcov

             [,1]         [,2]
[1,] 1.316277e+00 3.117570e-16
[2,] 3.117570e-16 1.039551e+02

標準誤差は分散共分散行列の平方根なので、

SEs <- sqrt(diag(vcov))

4.2.1 z値

それぞれの係数mを、それぞれの標準誤差SEsで割ると、帰無仮説（\(H_0\)）が\(0\)のz値が出て来ます。

zv <- m / SEs

両側検定を行えば、パッケージが表示するP値になります。

zp <- 2*(1 - pnorm(abs(zv)))
sprintf("z値%2.3f P値%2.3f", zv, zp)

[1] "z値10.489 P値0.000" "z値3.873 P値0.000"

4.2.2 t値

昔の統計解析パッケージでは、最尤法の推定量にもt値とそのP値が計算されていました。

tv <- m / SEs
tp <- 2*(1 - pt(abs(tv), n - length(tv)))
sprintf("t値%2.3f P値%2.3f",tv,tp)

[1] "t値10.489 P値0.000" "t値3.873 P値0.001"

自由度は(サンプルサイズ−パラメーターの数)になります。ただし、自由度の導出方法は確認できませんでした。

標準誤差の不偏推定量を使っていないので適切な計算方法では無いかも知れませんが、昔の統計解析パッケージのマニュアルを読む限りは特別に何か調整しているわけでは無さそうでした。

4.3 Wald検定

複数のパラメーターの値の合計を検定するときは、制約なしモデルの推定結果から検定等計量が計算できるWald検定の使い勝手が、t検定と似ていてよいです。

試しに\(\mu + \sigma^2 = 25\)を検定してみましょう。

b <- matrix(m, 2, 1) # 推定された係数を行列にする
C <- matrix(c(1, 1), 1, 2) # 制約を行列にする
R <- matrix(c(9 + 4^2), 1, 1) # 帰無仮説はゼロでなくてもよい
ws <- t(C %*% b - R) %*% solve(C %*% vcov %*% t(C)) %*% (C %*% b - R)
# χ二乗検定の自由度は制約の行列の行数
format(pchisq(abs(ws), 1, lower.tail = FALSE), digits=5)

     [,1]       
[1,] "0.0097372"

二次形式が\(\chi^2\)分布に従うことから、数理的に正当化されます。

複数の仮説、例えば\(\mu = 9\)で\(\sigma^2 = 16\)を検定することもできます。

b <- matrix(m, 2, 1) # 推定された係数を行列にする
C <- matrix(c(1, 0, 0, 1), 2, 2) # 制約を行列にする
R <- matrix(c(9, 4^2), 2, 1) # 帰無仮説はゼロでなくてもよい
# 仮説
ws <- t(C %*% b - R) %*% solve(C %*% vcov %*% t(C)) %*% (C %*% b - R)
format(pchisq(abs(ws), 2, lower.tail = FALSE), digits=5)

     [,1]       
[1,] "0.0021305"

なお、vcovはクラスター頑強分散共分散行列であることが求められるので、パネルデータのときは注意してください。

5 Rの関数で推定

上の例では制約条件がありませんでした。つまり重回帰分析の係数は無く、単に平均と分散を求めただけになります。重回帰分析の場合でも、対数尤度関数のy-muの部分がy-b0-b1*x1-b2*x2などと置き換わるだけなので、概念的には簡単です。ただし推定するパラメーター数と同じだけ一階条件の式が増え、ヤコビアンはその二乗の大きさになる事に注意してください。とてもでは無いですが、上の例のようにパッケージ無しで計算するのには無理があります。

Rにはnlmとoptimと言う関数で関数を最小化（もしくは最大化）するパラメーターを求めることができます。試しにnlmを使ってみましょう。

# 対数尤度関数を定義
llf <- function(p){
    # p[1]がμ，p[2]がσ^2
    s <- sum((y - p[1])^2)
    n <- length(y)
    -n*log(2*pi)/2 -n*log(p[2])/2 - s/p[2]/2
}
# 基本的に目的関数と初期値を与えれば計算してくれる
# グラディエントとヘッシアンは数値微分で内部的に計算
# nlm関数は最小化するパラメーターを探すので、対数尤度関数の符号を反転して目的関数にしている
r_nlm <- nlm(function(theta){ -llf(theta) }, c(11, 16), hessian = TRUE)
r_nlm

$minimum
[1] 97.70822

$estimate
[1] 12.03450 39.48839

$gradient
[1] 1.468968e-06 9.565458e-07

$hessian
              [,1]          [,2]
[1,]  7.597170e-01 -1.149791e-05
[2,] -1.149791e-05  9.615614e-03

$code
[1] 1

$iterations
[1] 13

solve(r_nlm$hessian)

            [,1]         [,2]
[1,] 1.316279571 1.573947e-03
[2,] 0.001573947 1.039975e+02

optim関数は、引数controlに渡す値で、目的関数を最大化するパラメーターを探す事ができます。

r_optim <- optim(c(11, 16), llf, hessian = TRUE, method = "L-BFGS-B", control = list(fnscale = -1))

もの凄く楽ですね。今までは何であったのかと言わないでください。なお、対数尤度関数やデータの都合によってはグラディエントをユーザー定義関数として与えた方がよいです。

5.1 尤度比（LR）検定

最尤法ではWald検定以外にも、尤度比検定とラグランジェ乗数検定がよく使われます。その中で、よく考えるとWilks定理が出てくるのでややこしい一方で、ぱっと見で直観的なのは尤度比検定です。帰無仮説\(\mu=10, \sigma=5\)を検定してみましょう。

# 尤度比検定の統計量; H0: μ=11, σ=6
lrs <- -2*(llf(c(10, 5^2)) - llf(r_optim$par))
# χ二乗分布からp値を求める; パラメーター2つを同時検定なので自由度は2
format(pchisq(lrs, 2, lower.tail = FALSE), digits=5)

[1] "0.013305"

上の例では不要でしたが、一般には制約無しのモデルと、制約付きのモデルをそれぞれ推定する必要があります。制約されていないパラメーターの値が分からないからです。制約に\(\mu=10\)を置いて、尤度比検定をかけてみましょう。

# 対数尤度関数をカリー化して、制約μ=10を置いて最適化
r_optim_restricted <- optim(c(11), function(p) llf(c(10, p)), hessian = TRUE, method = "L-BFGS-B", control = list(fnscale = -1))
# 尤度比検定の統計量; H0: μ=11, σ=6
lrs <- -2*(llf(c(10, r_optim_restricted$par)) - llf(r_optim$par))
# χ二乗分布からp値を求める; パラメーター1つを検定なので自由度は1
format(pchisq(lrs, 1, lower.tail = FALSE), digits=5)

[1] "0.083754"

lmやoptimに慣れれば尤度比検定をかけるのが楽です。

似たような検定が3種類あっても感があるわけですが、推定を2回する必要があるのが欠点

5.2 ラグランジェ乗数（LM）検定

スコア検定とも呼ばれるもので、直観的にはややこしそうですが、説明はそうでもないものです。制約付きモデルの推定と、グラディエントとヘッシアンの計算で検定統計量が計算できます。

# 数値微分を計算してくれるパッケージ
library("numDeriv")
# 制約付きモデルの推定結果から、制約付きパラメーターをつくる
# 今回の例は数学的にグラディエントとヘッシアンが導出してあるので、
# g <- matrix(c(eval(g1), eval(g2)), 2, 1)
# h <- matrix(c(eval(h11), eval(h21), eval(h12), eval(h22)), 2, 2)
# としても良いが、一般的には数値微分してしまうほうが手っ取り早いし、間違いが少ない
param_restricted <- c(10, r_optim_restricted$par)
# 制約なしモデルに、制約付きパラメーターを入れて、グラディエントをつくる
g <- matrix(grad(llf, param_restricted), 1)
# 制約なしモデルに、制約付きパラメーターを入れて、ヘッシアンをつくる
h <- hessian(llf, param_restricted)
# -hはフィッシャー情報行列
lmt <- g %*% solve(-h) %*% t(g)
# 制約されたパラメーターは1つなので、自由度1のχ二乗検定
format(pchisq(lmt, 1, lower.tail = FALSE), digits=5)

     [,1]      
[1,] "0.060894"

5.3 Wald検定，LM検定，LR検定の使い分け

この3つは本質的には同一の検定である事が示されているわけですが、Wald検定は制約無しモデルを推定すれば利用でき、LM検定は制約付きモデルを推定すれば利用でき、LR検定は2つのモデルをそれぞれ推定しないと利用できないです。計算量から言えばWald検定が優れているような気がしますが、現在の計算機の速度から言えば誤差ですし、コーディング量はLR検定が少ないです。Breusch-Pagan検定のような定番の手法で使っている場合でなければ、説明がしやすい手法を選びましょう。計量経済学を履修していれば、全部、試験範囲だと思いますが。

6 信頼区間

最尤法はヘッシアンから計算した標準誤差を元にしたワルド信頼区間と、尤度比検定の棄却域に境界が対応するプロファイル尤度信頼区間の二つが主に使われています。

6.1 ワルド信頼区間

漸近正規性を前提としている推定量です。

a <- 0.05 # 有意水準
lwrupr <- c(a/2, 1 - a/2)
estimated <- r_nlm$estimate
# Wald型信頼区間（左右対象）
confint.wald <- t(sapply(1:length(estimated), \(i) qnorm(lwrupr, mean = estimated[i], sd = SEs[i])))
rownames(confint.wald) <- c("μ", "σ²")
colnames(confint.wald) <- sprintf("%.1f%%", lwrupr*100)
print(confint.wald)

         2.5%    97.5%
μ   9.785854 14.28315
σ² 19.504913 59.47186

6.2 プロファイル尤度信頼区間

プロファイル尤度信頼区間は、計算に手間暇がかかるためか21世紀になってから一般化した信頼区間³で、今ではデフォルト出力になっていたりします。

尤度関数の複数の説明変数を、ひとつを残して陰関数として扱うことで、尤度関数を一変数関数にしたものがプロファイル尤度関数です。ある変数以外の変数は、ある変数にあわせて最適な値がとられるときの尤度が得られます。

計算には尤度比検定に用いる制約付き最尤法に使う関数が必要になります。グラディエント、ヘッシアン、制約付き対数尤度関数、制約付きグラディエントを用意しましょう。

# グラディエント
llfg <- function(p){
  n <- length(y)
  mu <- p[1]
  s2 <- p[2]
  c(sum(y-mu)/s2, -n/(2*s2) + sum((y-mu)^2)/(2*(s2^2)))
}
# ヘッシアン
hessian <- function(p){
  n <- length(y)
  mu <- p[1]
  s2 <- p[2]
  matrix(c(-n/s2, -sum(y - mu)/(s2^2), -sum(y - mu)/(s2^2), n/(2*(s2^2)) - sum( (y - mu)^2 )/(s2^3)), 2, 2)
}
pparam <- function(cp, i, v){
    # パラメーターのi番目をvとし、i番目以外はcpの値をcpの並びで使う
    p <- numeric(length(cp) + 1)
    if(1 < i) p[1:(i-1)] <- cp[1:(i-1)]
    p[i] <- v
    if(i < length(p)) p[(i + 1):length(p)] <- cp[i:length(cp)]
    p
}
pllf <- function(cp, i, v){
    # 制約をパラメーターに加えて対数尤度関数を呼ぶ
    llf(pparam(cp, i, v))
}
pllfg <- function(cp, i, v){
    # 制約をパラメーターに加えてグラディエントを計算し、制約部分以外に対応する要素を戻す
    llfg(pparam(cp, i, v))[-i]
}

まず、最尤法の結果から、その条件となるプロファイル尤度信頼区間の端点における対数尤度を求めます。

lwrupr_q <- qchisq(1 - a, df = 1, lower.tail = TRUE)
# 対数尤度がこの値の係数は、信頼区間の端点
lwrupr_llf <- -r_nlm$minimum - lwrupr_q/2 # r_nlm$minimumは対数尤度関数の符号が逆の値

プロファイル尤度信頼区間の端点にはもう一つ条件があり、信頼区間を計算するパラメーター以外のパラメーターのグラディエントはゼロです。

パラメーターを、対数尤度と信頼区間を計算する変数以外の変数のグラディエントに移すベクトル値関数を考え、これをニュートン・ラフソン法で推定することになります。信頼区間の上限と下限だけに解は二つあるのですが、最尤推定量のグラディエントはゼロで領域が分割されるため、初期値を最尤推定量の左側にするか右側にするかで、上限と下限をそれぞれ推定できます。

lsearch <- function(i, init.p){
    r_optim <- optim(rep(1, length(estimated) - 1), pllf, pllfg, i = i, v = init.p, method = "L-BFGS-B", control = list(fnscale = -1)) # methodをBFGSにすると収束しない
    b <- pparam(r_optim$par, i, init.p)
    for(j in 1:20){
        b0 <- b
        # グラディエントを一行変えて、ベクトル値関数の値にする
        f <- llfg(b)
        f[i] <- llf(b) - lwrupr_llf
        # ヘッシアンを一行変えて、ベクトル値関数のヤコビアンにする
        J <- hessian(b)
        J[i, ] <- llfg(b)
        b <- b - solve(J, f) # = solve(J) %*% f
        if(all(b - b0 < 1e-6)) break
    }
    b[i]
}
confint.lr <- matrix(NA, length(estimated), 2)
for(i in 1:length(estimated)){
    confint.lr[i, ] <- c(lsearch(i, estimated[i] - 2*SEs[i]), lsearch(i, estimated[i] + 2*SEs[i]))
}
rownames(confint.lr) <- c("μ", "σ²")
colnames(confint.lr) <- sprintf("%.1f%%", lwrupr*100)
(confint.lr)

         2.5%    97.5%
μ   9.711917 14.35710
σ² 24.760790 68.62141

指数分布族で制約が線形など典型的な場合はGLMを用いることで、収束しやすくなり、標準誤差が小さくなります。↩︎
Greene (2003) “Econometric Analysis,” 5th Edition↩︎
伊藤 (2019)に詳しい説明があります。↩︎

Rでべたべたと最尤法