今週も特にありません

進捗どうですか?

Rにおける最適化アルゴリズムをこれ一つで?

ほぼ使われているところを見ないoptimxパッケージを使ってみたメモ

データ分析を行っている方は、日々パラメータの推定などで最適化に励んでいると思います

Rで目的関数を自分で定義して最適化しようとする場合、optimやnlminb(一変数の場合はoptimize)等を用いて最適化を行います。 optimで最適化する場合、とりあえずデフォルトで実行されるNelder-Mead法を試して、この最適化が上手くいかなかった場合に、BFGS法などのその他の方法を試すという状況がよくあると思います

optimxパッケージは、Rで実行可能な最適化手法をまとめて一通り実行して、その結果を比較できるものです

尤度関数を定義して、そのパラメータを推定する例をやってみます

library(optimx)

loglik <- function(x, param) {
  -sum(dnbinom(x, size = param[1], mu = param[2], log = TRUE))
}

set.seed(71)
data <- rnbinom(500, size = 10, mu = 5)
initparam <- c(10, 5)
optimx(par = initparam, fn = loglik, x = data, control = list(all.methods = TRUE))

いくつかの警告文が出てきますが、最適化結果は以下のようになります

                   p1       p2         value fevals gevals niter convcode kkt1 kkt2 xtimes
BFGS         9.889487 5.038007  1.190725e+03     14      5    NA        0   NA   NA   0.01
CG           9.892441 5.037995  1.190725e+03    305    101    NA        1   NA   NA   0.27
Nelder-Mead  9.892127 5.038049  1.190725e+03     41     NA    NA        0   NA   NA   0.02
L-BFGS-B     9.889511 5.038000  1.190725e+03     13     13    NA        0   NA   NA   0.03
nlm          9.889480 5.037999  1.190725e+03     NA     NA    11        0   NA   NA   0.01
nlminb       9.889512 5.038000  1.190725e+03      8     18     5        0   NA   NA   0.00
spg          9.889768 5.038000  1.190725e+03     24     NA    18        0   NA   NA   0.13
ucminf       9.889507 5.037997  1.190725e+03      9      9    NA        0   NA   NA   0.02
Rcgmin             NA       NA 8.988466e+307     NA     NA    NA     9999   NA   NA   0.00
Rvmmin             NA       NA 8.988466e+307     NA     NA    NA     9999   NA   NA   0.00
newuoa       9.889511 5.038000  1.190725e+03     34     NA    NA        0   NA   NA   0.02
bobyqa       9.889509 5.038000  1.190725e+03     27     NA    NA        0   NA   NA   0.01
nmkb         9.890341 5.037991  1.190725e+03     69     NA    NA        0   NA   NA   0.03
hjkb        10.000000 5.000000  1.190775e+03      1     NA     0     9999   NA   NA   0.02

関数自体の使い方は普通のoptimとほとんど変わりません。 all.methods=TRUEを指定することですべての手法を実行します。 p1、p2が最適化したパラメータの値、valueがその時の目的関数の値で、 convcodeが0以外のものは最適化が正常に終了していないことを表しています

今回の場合では、結果が大きく変わることはないですが、 より複雑な場合だとそれぞれのアルゴリズムを一覧で比較できるところは使いどころかもしれません

より詳しい説明はJournal of statistical softwareにもあります。 CRAN Task View: Optimization and Mathematical Programmingを見ると、知らないパッケージがたくさんあるなー、まだまだRはいろいろとできるのだなーと感じます