関数の最大・最小化の編集

//間瀬茂
COLOR(red){SIZE(30){汎用最適化関数 optim() の使用法}}

-統計学では、最尤推定や最小自乗法等、ある目的関数の最適化を行なう必要が
しばしば起こる。COLOR(red){optim()} は R の汎用最適化関数で、複数の代表的手法をオプションで選べる便利な道具である。
-一変数関数専用の最適化関数 COLOR(red){optimize} もある。
-線形不等式制約の場合は [[線形不等式制約付きの最適化関数 constrOptim]] 参照。
-同種の関数として [[汎用非線形最小化関数 nlm]] もある。
-一変数関数の根を求める専用関数 COLOR(red){uniroot} もある。
-非線形回帰専用の関数 COLOR(red){nls} は独自の最適化アルゴリズムを使用する。

#contents

*SIZE(20){optim()} の使用法 [#z6cdf87c]

**命令書式 [#sc95e06a]

optim(par, fn, gr = NULL,
       method = c("Nelder-Mead", "BFGS", "CG", "L-BFGS-B", "SANN"),
       lower = -Inf, upper = Inf,
       control = list(), hessian = FALSE, ...)
 
COLOR(red){optim()} は再帰的に使うことができる(目的関数自身が内部で COLOR(red){optim}  を使っていても良い）。引数 COLOR(red){par} はスカラーでも、ベクトル(多変数関数)でも良い。

**引数の意味 [#w4e72c3a]

*** par [#yc3cb4c4]
目的関数のベクトル引数に対する初期値。この選択は一般に試行錯誤によらざるを得ない。まずい選択をすると。収束しないか、局所的最的値しか得られない。

*** fn [#o8af719e]
最小化(既定)もしくは最大化すべき目的変数。最初の(ベクトル)引数に関して最適化が実行される。返り値はスカラーである必要。COLOR(red){fn} は特定のパラメータ値で関評価不能なら COLOR(red){NA} や  COLOR(red){Inf} 値を返しても良いが、初期パラメータ値では有限の値を返さなければならない(COLOR(magenta){"L-BFGS-B" 法}では常に有限の値を返す必要)。

*** gr [#bc97c2ae]
"BFGS", "CG", "L-BFGS-B" 法に対するグラディエント(一階(偏)微分)関数。もし COLOR(red){NULL} を指定(既定)すると、グラディエント関数は微分商を用いた数値微分で計算される。もし、比較的簡単に定義できるなら明示的に与えた方が良い。"SANN" 法に対しては、これは新しい候補パラメータを計算する関数と解釈され、もし COLOR(red){NULL} なら既定の Gaussian Markov kernel が使用される

*** method [#sf1e857d]
使用する最適化手法を指示する。既定手法は COLOR(magenta){"Nelder-Mead" 法}

*** lower,  upper [#h6b41a03]
COLOR(magenta){"L-BFGS-B" 法}に対する変数の下限、上限を与える。既定値は COLOR(red){lower = -Inf}, COLOR(red){upper = Inf}

*** control [#a8dbd5e2]
制御パラメータのリスト。

*** hessian [#z8ae39f1]
論理値。数値的に計算したヘッシアン(二階偏微分関数行列の値)を返すべきかどうかを指示。

*** ... [#r9bd65a5]
COLOR(red){fn} と COLOR(red){gr} に引き渡される追加の引数(群)

**使える最適化手法 [#j7aa73f2]

*** "Nelder-Mead" 法 [#a0920368]
COLOR(magenta){Nelder-Mead 法}。関数値だけを用い、頑健(例えば初期値の選択に敏感でない)であるが、相対的に遅い。微分できない関数に対してもそれなりに使える

*** "BFGS" 法 [#h87459ad]
COLOR(magenta){準ニュートン法}。COLOR(magenta){variable metric 法}とも呼ばれる。関数値とグラディエント関数を最適化関数の曲面を近似するのに使う。Broyden, Fletcher, Goldfarb, Shanno の4人が同時期に提案。

*** "CG" 法 [#k48c0acb]
Fletcher and Reeves によるCOLOR(magenta){共役勾配法}。Polak-Ribiere と Beale-Sorenson による更新法も選択できる共役勾配法は準ニュートン法に比べると破綻しやすいが、メモリー使用量が少ないため、より大規模の最適化に使える可能性

*** "L-BFGS-B" 法 [#c303cc6e]
各変数が上限・下限によるCOLOR(magenta){制約条件を許す準ニュートン法}の変種。初期値はこの制約条件を満たす必要

*** "SANN" 法 [#z688fed4]
確率的手法であるいわゆるCOLOR(magenta){シミュレーテッド・アニーリング法}。関数値だけを用い遅い。微分できない関数にも使える。COLOR(magenta){メトロポリス法}を用いる。COLOR(magenta){組合せ的最適化問題}(巨大な有限集合内から最適値を探す)にも使える。一般的手法ではないが、非常にラフな関数に対しても良い値を得ることができる可能性がある

** 制御パラメータリスト control の意味 [#uaacef6e]

*** trace [#j62e2900]
非負整数。もし正なら途中結果が表示される。値が大きいほどより詳しい情報

*** fnscal [#gf14ec83]
最適化の途中で、関数 COLOR(red){fn} と COLOR(red){gr} に適用される比例定数。最適化は実際は変数 COLOR(red){fn(par)/fnscale} に対して行われる。もし負ならば、最適化は実は最大化となる

*** parscale [#h3018f50]
パラメータに対する比例定数(ベクトル)。最適化の際のパラメータは実際は COLOR(red){par/parscale} とされる(パラメータの変化を誇張・抑制できる)。

* 注意 [#e28ecc95]

** 最大化するには？ [#ecf7726d]
SIZE(20){optim()} は既定では最小化を行なう。最大化をするには制御リスト COLOR(red){control} の COLOR(red){fnscale} に負の値を与える(例えば COLOR(red){fnscale=-1})。本来の目的関数の代わりに、値にマイナスを掛けた目的関数を用意するのもありです。

** 非線形連立方程式を解くには？ [#h9b8656f]
*** optim関数 [#hb12e05f]
非線形連立方程式 COLOR(red){F1(x)=F2(x)=...=Fn(x)=0} を解くには、例えば目的関数を COLOR(red){F(x)=F1(x)^2+F2(x)^+...+Fn(x)^2} と置いて COLOR(red){F(x)} の最小化を行なう。
*** RMinpackパッケージ [#s8a530f2]
http://www.hppi.troitsk.ru/Kondrin/ にあるRMinpack(CRANにはまだないパッケージ)のFSolve関数を用いる。

*** nleqslvパッケージ [#v7f2c799]
[[連立方程式を解く]]を参照

** 非線形回帰をするには? [#gde33b2f]
データ COLOR(red){d[i]} にパラメータベクトル COLOR(red){t} を持つ非線形関数 COLOR(red){f[i]=F(i, t)} を当てはめるには、目的関数を 
COLOR(red){SS(t)=sum((d-f)^2)} と置いて COLOR(red){t} について最小化すれば良い。

** 変数に制約がある時は? [#i1f12672]
-変数に制約がある時は最適化は一層困難になります。矩形型の制約なら "L-BFGS-B" 法が使えるかもしれません。
-しばしば使われる方法は、制約条件を表す関数 COLOR(red){g(x)} (例えばこの値が大きい COLOR(red){x} ほど望ましくない) を用い、目的関数を COLOR(red){f(x)} から COLOR(red){f(x)+g(x)} に変えることです。真の制約条件関数 COLOR(red){g(x)} は無限大の値を取る(不可能を意味する)のが普通でしょうが、そうすると例えば微分できなくなりますから、無限大の代わりに相当大きな値をとるように便宜的に定義します。
-線形不等式制約の場合は [[線形不等式制約付きの最適化関数 constrOptim]] 参照。

**本当に最適値が求まったの? [#ybf88577]
-これは意外に難しい点です。如何なる数値的最適化手法も、COLOR(magenta){原則として局所的最適値しか見つけられない}と心得るべきです。
-さらに計算に於けるCOLOR(magenta){数値誤差の問題}も無視できません。
-普通COLOR(magneta){最適値があると信じ}て計算するのでしょうが、本当でしょうか？
-また、本当に真の最適値が必要なのでしょうか(良くある例では、COLOR(magenta){目的関数がかなり小さくなればそれで良し})。
-まず収束したかどうかを結果の収束判定コード COLOR(red){convergence} が COLOR(magenta){0} かどうかでチェックします。もし収束していてもすぐ信用してはいけません。初期値を変えて、そして使う手法を変え、何度も計算し、結果を比較すべきです(普通微妙に、時に信じられないほど異なります)。
-収束判定用の基準値 COLOR(red){abstol}, COLOR(red){reltol} 等を変えると結果も変わるかもしれません。
-もし、収束しないままに打ち切られていたら、たとえ、尤もらしい COLOR(red){par} 値が与えられて信用できません。
-微分可能な関数ならば、(制約付き最適化は例外かもしれない)理論的には最適値ではグラディエントベクトルはゼロベクトルになるはずですが、普通は微妙に異なります。その異なり方を子細に見るのも役に立ちます。
-いろいろ書くとやる気を無くしそうですが、心配しないで下さい。ほとんどの場合、皆さん出てきた結果が真の最適値と信じている、そしてそれはそれで何とかなっている(らしい)ことが事実なのですから。
-統計学ではランダムなデータを扱いますから、誤差が大きければ(そしてデータ数が少なければ)、たとえ真の値がわかっていても、最適化の結果は大きくずれることがあっても当然です。

** 初期値パラメータの選び方は？ [#kbc1995a]
***神頼み、人頼み [#d73dac7a]
-基本的に試行錯誤を繰り返して選ぶより仕方がないと心得る
-問題の性質からありそうな範囲の大雑把な見当をつけるべきである
-教科書には書いてないお勧めの方法として、COLOR(magenta){身近の偉い人、経験者}に相談するのも良い方法である。

***システマティックな探索 [#ae232611]
-パラメータの範囲を適当に等分して、各升目毎に、目的関数の値を網羅的に計算する。その中から最適値に近そうなパラメータ値を選ぶ。
-1変数ならば COLOR(red){plot()} 関数でグラフを描いてみる。
-2変数ならば、計算結果を行列に表現した上で、 COLOR(red){image()} 関数や COLOR(red){contour()} 関数で視覚化してみるとわかりやすいかもしれない。

*** 多段階最適化 [#n455d54b]
-COLOR(magenta){"Nelder-Mead" 法} や COLOR(magenta){"SANN" 法} は相対的に時間がかかるが、初期値の選択に比較的敏感でない。これらをまず使い、初期値の候補を求める。
-次に、それを初期値として、改めて他の切れ味の鋭い(ただし刃こぼれしやすい)方法を使う。COLOR(red){optim(optim(x, fn)$par, fn, method="BFGS") }といった使い方もできます。
-微分できないような関数ではCOLOR(magenta){"Nelder-Mead" 法} やCOLOR(magenta){"SANN" 法} しか選択肢は無い。

*** それでも駄目なら (1) [#ta9ecccd]

-問題の設定を変えてみる。問題を同値なより単純な問題に整理する。目的関数を変形する(対数を取る、分母を払う)。
-できれば、変数の数を減らした問題に置き換える。
-変数を置き換えて制約条件を除く。例えば、確率パラメータ COLOR(red){p} は COLOR(red){p=1/(1+e^x)}) と置き換えれば制約がなくなる。正のパラメータ COLOR(red){x} は COLOR(red){x=e^y} と置けば制約がなくなる。
-少しは頭を使い、実は紙と鉛筆で解ける問題ではないか反省してみる。
-繰り返し回数を増やしてみる。
-COLOR(red){parscale} を大きくし、一回毎のパラメータの変化を小さくしてみる。
-COLOR(red){fnscale} を小さな(正)数にし、関数値の凹凸を誇張する。-COLOR(red){trace} オプションを付けて途中経過を眺めると、問題点がはっきりすることがある。

*** それでも駄目なら (2) [#h8e5e6fb]
目的関数をある手法で最適化した際の、COLOR(red){返り値の適当な関数}を目的関数とし、関数値のみを使用する手法である COLOR(magenta){"Nelder-Mead" 法} で再帰的に最適化を行なう。例えば、非常に平坦な関数で、関数値だけではなかなか最適値が決まらないとする。最適値ではグラディエントベクトルは理論的にゼロベクトルになるべきであるから、その自乗和を目的関数に取ることが考えられる。
-optim(par, fn=optim(par, fn=SS)$value) # ただし非常に時間がかかる。

*** それでも駄目なら (3) [#s265c525]
人生には最適化問題よりもっと大事なことがあると達観し、さっさとあきらめた方が良いかも。すべての最適化問題が解けるとは限りません。最尤推定が駄目なら、もっと簡単なモーメント型推定量を考える、等。もちろん R の COLOR(red){optim()} よりも優れたソフトがある可能性はいつでもあります。
[[汎用非線形最小化関数 nlm]] も試してみて下さい。しばしば良い結果を得るという報告があります。

*実例 [#df5368fe]

** 実例 (1) 負の二項分布のパラメータの最尤推定 [#s8619974]

既定の Nelder-Mead 法で対数尤度を最大化(オプション fnscale=-1  を忘れず加える)

d <- c(1612,164,71,47,28,17,12,12,5,7,6,3,3,13) # データ(最後は個数13以上)
 NB <- function (s, p) { # パラメータ (size, prob)=(s, p) の負の二項分布の確率関数
   P <- dnbinom(0:12, size=s, prob=p, log=T) # 個数 0,1,...,12 の確率の対数
   P[14] <- pnbinom(12, size=s, prob=p, lower.tail=F, log.p=T) # 個数13以上の確率の対数
   return(P)
 }
 LL <- function (x) return(sum(NB(x[1],x[2])*d)) # 目的関数(対数尤度を計算する関数)
 optim(c(0.2, 0.5), LL, control=list(fnscale=-1)) # 初期パラメータ (0.2, 0.5) で最大化実行
 $par
 [1] 0.1155133 0.1553361  # (size, prob) パラメータの最尤推定値
 $value
 [1] -1705.324            # その時の対数尤度値
 $counts
 function gradient
       63       NA        # Nelder-Mead 法繰り返し数63回、gradient 評価は無し
 $convergence
 [1] 0                    # 収束判定コード 0 (無事収束)
 $message
 NULL                     # その他のメッセージ無し
 > sum(d)*exp(NB(res$par[1],res$par[2])) # 推定度数の計算
  [1] 1612.914091  157.371829   74.140527   44.160512   29.052777   20.198801
 [7]   14.546131   10.734116    8.064296    6.142199    4.729214    3.673399
 [13]    2.874090   11.398017

**実例 (2) COLOR(red){example(optim)} から [#a2520f03]

> fr <- function(x) { # 目的関数 (Rosenbrock の Banana 関数)
    x1 <- x[1];  x2 <- x[2]
    100 * (x2 - x1 * x1)^2 + (1 - x1)^2
 }
 grr <- function(x) { # そのグラディエント関数 (返り値がベクトルに注意)
    x1 <- x[1]; x2 <- x[2]
    c(-400 * x1 * (x2 - x1 * x1) - 2 * (1 - x1), 200 * (x2 -x1 * x1))
 }
 # 以下各種手法による最適パラメータと最適値を紹介
 # 多数決により、最適値は 1,1 らしい(!?)
 # でも最適値は結構違う(ように見える)
 > optim(c(-1.2, 1),  fr)  # 既定の "Nelder-Mead" 法
 $par
 [1] 1.000260 1.000506
 $value
 [1] 8.825241e-08
 > optim(c(-1.2, 1),  fr,  grr,  method = "BFGS")  # "BFGS" 法
 $par
 [1] 1 1
 $value
 [1] 9.594955e-18
 > optim(c(-1.2, 1), fr, NULL, method = "BFGS", hessian = TRUE)
 $par
 [1] 0.9998044 0.9996084
 $value
 [1] 3.827383e-08
 > optim(c(-1.2, 1), fr, grr, method = "CG")
 $par
 [1] 1 1
 $value
 [1] 4.477649e-17
 > optim(c(-1.2, 1), fr, grr, method = "CG", control = list(type = 2))
 $par
 [1] 1 1
 $value
 [1] 7.12628e-18
 > optim(c(-1.2, 1), fr, grr, method = "L-BFGS-B")
 $par
 [1] 0.9999997 0.9999995
 $value
 [1] 2.267630e-13

**実例 (3) COLOR(red){example(optim)} から [#s8804044]

> fw <- function(x) 10 * sin(0.3 * x) * sin(1.3 * x^2) +  # あるワイルドな目的関数
                              1e-05 * x^4 + 0.2 * x + 80
 > plot(fw, -50, 50, n = 1000)  # グラフを書いてみると様子がよくわかる
 > res <- optim(50, fw, method = "SANN", control = list(maxit = 20000,
                        temp = 20, parscale = 20))  # "SANN" 法で第一段階最適化
 > res               # その結果
 $par
 [1] -15.81488  #  # 真のパラメータは約 -15.81515
 $value
 [1] 67.46834
 > optim(res$par, fw, method = "BFGS")  # "BFGS" 法で第二段階最適化
 $par
 [1] -15.81515   # "真のパラメータ" が求まった！
 $value
 [1] 67.46773    # 最適値も確かに前より小さい

**実例 (4) COLOR(red){非線形回帰 (日本人の名字100位の頻度への Zipf 分布の当てはめ)} [#y06abdd6]
 # 日本人の代表的苗字100位(電子電話帳を集計したデータ。単位世帯数)
 d <- c(456430,403506,335288,314770,256706,255876,254662,249509,241651,203101,
       197460,193503,169617,152065,149006,143552,137475,137160,129673,123953,
       114802,110430,108369,108345,105778,102647, 97704, 95699, 93207, 91298,
        90925, 89856, 89818, 87815, 86992, 86695, 86234, 78849, 78178, 76233,
        75826, 75264, 74510, 74352, 73185, 72640, 72569, 71443, 71102, 70889,
        70082, 69904, 68661, 67852, 67571, 65830, 64234, 64119, 63180, 60606,
        59967, 58141, 57568, 57037, 56651, 56538, 56324, 55793, 55208, 54878,
        53284, 52858, 50891, 50499, 50349, 50180, 49474, 49397, 49026, 48747,
        48724, 48386, 48329, 47744, 47094, 46923, 46858, 46621, 46098, 46004,
        45682, 45164, 44731, 44650, 44641, 44222, 43908, 43763, 43669, 43205)
 d <- d/sum(d)  # 頻度
 zipf <- function (x) { # Zipf 分布 p(i) = c/(i^x), i=1,2,...,100
   s <- 1/(1:100)^x
   return(s/sum(s))}   # 最後に正規化
 SS <- function (x)  sum((d-zipf(x[1]))^2)   # 目的関数 （誤差の自乗和）
 res <- optim(0.7,  SS,  control=list(trace=TRUE, parscale=10))  # Nelder-Mead 法で当てはめ
 > res$par
 [1] 1.030962  # 非線形最小自乗推定値
 > 1-SS(res$par)/sum(d^2)
 [1] 0.9336638       # 決定係数 (結構高い)

タイムスタンプを変更しない