対数尤度関数の最適化

汎用最適化関数 optim

R の optimize 関数が変数を 1 つだけもつ目的関数の最適化に用いるのに対し、変数が 2 つ以上のときも利用できるのが optim 関数である。

正規分布のパラメーター推定

正規分布のパラメーターとして平均と分散の 2 つある。optim 関数を利用して、最尤推定法によりパラメーターを推定する例を示す。まず、平均 10、分散 5 の正規分布に従う乱数を 100 個生成し、変数 y に格納する。

y <- rnorm(100, mean = 10, sd = 4)

次に正規分布のパラメーター推定時に用いる対数尤度関数を作成する。パラメーターを x とし、x[1] を平均、x[2] を分散とする。

loglikelihood <- function(x, y) {
  - sum(0.5 * (y - x[1]) ^ 2 / x[2] + 0.5 * log(x[2]))
}

optim 関数を利用して対数尤度関数を最大にする x[1]x[2] を求める。推定結果をみると、x[1] が 9.5、x[2] が 15.4 と乱数生成条件に近いことが確認できる。optim 関数はデフォルトで最小化を行うため、control 引数に fnscale = -1 を代入して、最大化を行うようにする。

optim(c(10, 2), loglikelihood, y = y, control = list(fnscale = -1))
## $par
## [1]  9.533245 15.429277
## $value
## [1] 186.7908
## $counts
## function gradient 
##       55       NA 
## $convergence
## [1] 0
## $message
## NULL

負の二項分布のパラメーター推定

負の二項分布のパラメーターを最尤推定法により推定する例を示す。対数尤度関数の最大化は optim 関数を利用する。尤度関数と確率分布関数が同じであることを利用して、ここでは尤度関数を自作せずに R で用意された負の二項分布の確率分布関数 dnbinom を利用する。

y <- rnbinom(100, mu = 1000, size = 200)

loglikelihood <- function(x, y){
  sum(log(dnbinom(y, mu = x[1], size = x[2])))
}

optim 関数を利用してパラメーター推定を行う(対数尤度関数の最大化を行う)。optim 関数はデフォルトで最小化を行うため、control 引数に fnscale = -1 を代入して、最大化を行うようにする。

optim(c(900, 100), loglikelihood, y = y, control = list(fnscale = -1))
## $par
## [1] 995.1058 193.4067
## $value
## [1] 577.6733
## $counts
## function gradient 
##       63       NA 
## $convergence
## [1] 0
## $message
## NULL

最適化アルゴリズム

R の optim 関数を利用して最適化を行う際に、手法の指定として Nelder-Mead、BFGS、CG、L-BFGS-B、SANN を指定することができる。これらの手法は以下のような特徴を持つ。

Nelder-Mead 目的関数の変数が n 個の場合に、n + 1 個の頂点を持つ多面体を作成し、n + 1 面体の頂点での目的関数を評価することによって最適化を行う。評価回数が多いために計算に時間がかかる。しかし、導関数がわからなくても利用できるというメリットがある。
BFGS 準ニュートン法とも呼ばれ、関数の極大極小解を探索する方法である。
CG 共役勾配法という。連立一次方程式を解くための方法であるが、目的関数の最小化などの最適化問題を解くために用いることもできる。
L-BFGS-B
SANN

References

  1. Example of MLE Computations, using R. PDF
  2. Lecture 7: Negative binomial model. HTML 2010.