計量分析で最初に習うOLSは性質がよく計算しやすい推定ですが、非線形モデルの推定に使えません。非線形関数を推定する場合は、その非線形関数を推定するパッケージを用いない限りは、最尤法が手っ取り早い選択になります1。最尤法は、一般化線形回帰モデル（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元の連立方程式をベクトル値関数\(f(\mu, \sigma^2)\)にまとめます。
\[ \begin{aligned} f(\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} f(\mu, \sigma^2) = f(\mu_0, \sigma_0^2) + \frac{\partial f}{\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 f}{\partial(\mu_0, \sigma_0^2)}^{-1} f(\mu_0, \sigma_0^2) \end{aligned} \]
\(f\)が凸関数であれば、右辺の値は\(f(\mu_0, \sigma_0^2)\)よりも\(f(\mu, \sigma^2)\)に近づくので、右辺の計算を繰り返し行なえば解になる不動点に辿り着けます。
3.2 実際に計算するコード
上述の説明をだいたい忠実に実装したコードは以下になります。
# サンプル
set.seed(31)
y <- rnorm(30, mean=12, sd=7)
n <- length(y)
# 連立方程式を設定
f1 <- expression(-sum(y-mu)/s2)
f2 <- expression(-n/(2*s2) + sum((y-mu)^2)/(2*(s2^2)))
# muとs2で一階微分を作っておく
g11 <- expression(n/s2)
g12 <- expression(sum(y - mu)/(s2^2))
g21 <- expression(sum(mu - y)/(s2^2))
g22 <- 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)で評価したf1とf2の値を2x1行列にする
f <- matrix(c(eval(f1), eval(f2)), 2, 1)
# ヤコビアンに(mu, s2)を代入した行列を作る
j <- matrix(c(eval(g11), eval(g21), eval(g12), eval(g22)), 2, 2)
# 行列mから、ヤコビアンの逆行列と評価式を乗じたものを引く
m <- m - solve(j)%*%f
# 行列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 統計学的な意味
パラメーターの値の次は、標準誤差を求ます。 まず、2つのパラメーター\(\mu, \sigma\)をベクトル\(\theta\)に置きなおします。
\[ \theta = \begin{bmatrix} \mu \\ \sigma^2 \end{bmatrix} \]
上の例のヤコビアン\(\partial f/\partial \theta\)は、対数尤度関数を2階偏微分した行列なので、対数尤度関数のヘッシアン\(H(\theta)\)になるのに注意して、真のパラメーター\(\theta_p\)の周りで、テイラーの定理を用います。2階微分の引数\(\bar{\theta}\)は、平均値の定理から存在することが分かる、\(\hat{\theta}\)と\(\theta_p\)の間にある値です。
\[ f(\hat{\theta}) = f(\theta_p) + H(\bar{\theta})(\hat{\theta} - \theta_p) = 0 \]
\[ \hat{\theta} - \theta_p = - H(\bar{\theta})^{-1} f(\theta_p) \]
両辺に\(\sqrt{n}\)を乗じます。
\[ \sqrt{n} (\hat{\theta} - \theta_p) = - H(\bar{\theta})^{-1} \sqrt{n} f(\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} f(\theta_p)\)に分布収束することが分かります。
\[ \sqrt{n} (\hat{\theta} - \theta_p) \xrightarrow{d} - H(\theta_p)^{-1} \sqrt{n} f(\theta_p) \]
\(f(\theta_p)\)に\(n\)で割り、\(H(\theta_p)^{-1}\)に\(n\)を乗じます。\(\bar{f}(\theta_p)=f(\theta_p)/n\)に注意すると、
\[ \sqrt{n} (\hat{\theta} - \theta_p) \xrightarrow{d} - \big[ \frac{1}{n} H(\theta_p) \big]^{-1} \sqrt{n} \bar{f}(\theta_p) \]
と書けます。
\(\sqrt{n}\bar{f}(\theta_p)\)は乱数の平均値になるので、Lindberg–Levy中心極限定理を用いることができ、
\[ \sqrt{n} \bar{f}(\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{f}(\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(abs(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.0088456"複数の仮説、例えば\(\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.42987"なお、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.467787e-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] 13solve(r_nlm$hessian) [,1] [,2]
[1,] 1.316279571 1.573947e-03
[2,] 0.001573947 1.039975e+02optim関数は、引数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(f1), eval(f2)), 2, 1)
# h <- matrix(c(eval(g11), eval(g21), eval(g12), eval(g22)), 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検定は制約無しモデルを推定すれば利用でき、LR検定は制約付きモデルを推定すれば利用でき、LR検定は2つのモデルをそれぞれ推定しないと利用できないです。計算量から言えばWald検定が優れているような気がしますが、現在の計算機の速度から言えば誤差ですし、コーディング量はLR検定が少ないです。Breusch-Pagan検定のような定番の手法で使っている場合でなければ、説明がしやすい手法を選びましょう。計量経済学を履修していれば、全部、試験範囲だと思いますが。
指数分布族で制約が線形など典型的な場合はGLMを用いることで、収束しやすくなり、標準誤差が小さくなります。↩︎