RでべたべたとBHHH法 - 餡子付゛録゛

Berndt-Hall-Hall-Hausmanアルゴリズムは、化石的な最適化アルゴリズムで、実証分析を行っている経済学徒でお世話になったことがある人は、そこそこいると思います。TSPと呼ばれる統計解析パッケージ*1では標準的に使えたのですが、RではmaxLikパッケージ*2などを入れる必要があります。利用困難と言うわけではないですが、人気ではなさそうです。マルコフ連鎖モンテカルロ法が当たり前の時代、古の技法になりつつあると言えるでしょう。しかし、ちょっとした都合でBHHH法の計算手順を確認したいと思います。

1. BHHH法の概要

目的関数として対数尤度関数 $\mathcal{l}(y; \beta)$ が与えられたときのことを考えましょう。 $y$ は被説明変数のベクトル、 $\beta$ はパラメーターです*3。もっとも基本的なものはニュートン・ラフソン法になります。真の係数を $\beta$ 、推定量を $b^{(m)}$ （ $m$ は反復回数で指数ではない）とすると、

$b^{(m)} = b^{(m-1)} - \alpha^{(m-1)} \frac{\partial \mathcal{l}^{(m-1)}}{\partial{b^{(m-1)}}\partial{b^{(m-1)}}} \frac{\partial{\mathcal{l}^{(m-1)}}}{\partial{b^{(m-1)}}}$

になります。ここで、二階微分の計算があると上手く計算されないことが多々あるので、一階微分しかないフィッシャー情報行列の符号を反転させたもので代用しようと言うのが、BHHHアルゴリズムです。原理上、最尤法の計算ぐらいにしか使えません。

2. データセット

サンプルサイズ10のデータセットです。

X <- t(matrix(c(1, 5.75095, 7.13628, 1, 8.12624, 3.43466, 1, 0.470068, 8.32101, 1, 6.97381, 8.56067, 1, 3.01102, 1.86405, 1, 2.73131, 9.86068, 1, 0.639832, 7.30122, 1, 2.02848, 3.03172, 1, 2.27631, 0.179666, 1, 2.59997, 7.71923), 3))
y <- t(matrix(c(-8.67053, 6.59171, -21.7837, -10.9077, 1.78061, -23.3735, -21.2741, -3.17336, 4.87379, -16.6442), 1))

3. 対数尤度関数とグラディエント

正規分布を仮定した線形回帰を最尤法で求めます。対数尤度関数を一階微分したグラディエントは、フィッシャー情報行列の計算のために、説明変数と従属変数を一つづつ取れるようにします。

# 対数尤度関数
objf <- function(p){
  theta <- p[1]
  len <- length(p)
  b <- matrix(p[2:len], len-1, 1)
  s <- sum((y - X %*% b)^2)
  n <- length(y)
  -n*log(theta^2)/2 - s/theta^2/2
}

# グラディエント
objg <- function(p, y=get("y"), X=get("X")){
  n <- length(y)
  theta <- p[1]
  len <- length(p)
  b <- matrix(p[2:len], len - 1, 1)
  res <- y - X %*% b
  dtheta <- -1*n/theta + sum(res^2)/theta^3
  dbeta <- numeric(length(b))
  for(i in 1:length(b)){
    dbeta[i] <- sum(X[, i]*2*res)/theta^2/2
  }
  c(dtheta, dbeta)
}

4. フィッシャー情報行列の計算

フィッシャー情報行列の計算にヘッシアンが要りそうな感じがする*4のですが、 $i$ 行 $j$ 列が $\frac{1}{N} \sum_{n=1}^N \frac{ \partial \mathcal{l}(y_n; \beta) }{ \partial \beta_i } \frac{ \partial \mathcal{l}(y_n; \beta) }{ \partial \beta_j }$ となる分散共分散行列なので素朴に計算できます。これを観測数で割ったもの*5に-1を乗じて、ヘッシアンの代わりにします。

FIM <- function(p){
  m <- matrix(0, length(y), length(p))
  for(i in 1:length(y)){
    m[i, ] <- objg(p, y[i], matrix(X[i, ], 1))
  }
  (t(m) %*% m)/nrow(m)
}

なお、サンプルサイズを増やすと各セルの値が小さくなりすぎて「システムは数値的に特異です: 条件数の逆数 = ･･･」とエラーが出たりするので、実用的に使うには、もっと改良が要ります。

5. ステップ幅の計算

ステップ幅 $\alpha$ の具体的な計算方法については原論文やアルゴリズム紹介には言及が無いので、「Rで黄金分割探索」で定義した関数を流用します。BFGS法の実装でも、ステップ幅の計算はまちまちのようで、プログラマに任されているのでしょう。

6. BHHH法による推定

ニュートン・ラフソン法とほとんど一緒です。

init <- c(1, 1, 1, 1) # 初期値
param <- init # 更新していくパラメーターの変数
pv <- objf(param) # 初期値の目的関数の値を保存

for(i in 1:300){ # 収束が遅いので繰り返し回数上限は長めに
  g <- matrix(objg(param, y, X), length(param))
  R <- -1*( FIM(param) / nrow(X) )
#
# 以下に置き換えるとニュートン・ラフソン法になる
#  library("numDeriv")
#  R <- jacobian(function(p){ objg(p, y, X) }, param)
#
  # ステップ幅を計算
  a <- gss(function(a){
    # solve(R) %*% g = solve(R, g)
    objf(param - a * solve(R, g))
  }, min=0, max=1, maximize=TRUE) # データセットを変えてみて、うまく収束しなくなった場合はmin=-1にすると･･･

  # パラメーターを更新
  param <- param - a * solve(R, g)

  # 新たなパラメーターで目的関数の値を計算
  v <- objf(param)

  # 改善幅が小さければ終了
  if(abs(pv - v) < 1e-12){
    print(sprintf("pv - v = %f", pv - v))
    break;
  }

  # 計算の途中経過をばらばらと表示
  print(sprintf("i:%d a=%f", i, a))
  print(round(c(param), 6))

  # 現在のパラメーターの目的関数の値を保存
  pv <- v
}

# 推定結果をリストにまとめる
r_bhhh <- list(
  param = param,
  g = matrix(objg(param, y, X), length(param)),
  R = -1*( FIM(param)/ nrow(X) ) ,
  vcov = -1*solve( FIM(param))
)

7. 推定結果の比較

係数はほぼ同じモノになります。

r_nlm <- nlm(function(p){ -1*objf(p) }, init, hessian=TRUE)
r_ols <- lm(y ~ 0 + X)
r <- matrix(c(r_bhhh$param, r_nlm$estimate, NA, coef(r_ols)), ncol(X) + 1)
colnames(r) <- c("BHHH", "nlm/DSM", "lm/OLS")
r