![[RjpWiki]](image/../Rlogo.png "[RjpWiki]")
汎用最適化関数 optim() の使用法
optim(par, fn, gr = NULL,
method = c("Nelder-Mead", "BFGS", "CG", "L-BFGS-B", "SANN"),
lower = -Inf, upper = Inf,
control = list(), hessian = FALSE, ...)
optim() は再帰的に使うことができる(目的関数自身が内部で optim を使っていても良い）。引数 par はスカラーでも、ベクトル(多変数関数)でも良い。
目的関数のベクトル引数に対する初期値。この選択は一般に試行錯誤によらざるを得ない。まずい選択をすると。収束しないか、局所的最的値しか得られない。
最小化(既定)もしくは最大化すべき目的変数。最初の(ベクトル)引数に関して最適化が実行される。返り値はスカラーである必要。fn は特定のパラメータ値で関評価不能なら NA や Inf 値を返しても良いが、初期パラメータ値では有限の値を返さなければならない("L-BFGS-B" 法では常に有限の値を返す必要)。
"BFGS", "CG", "L-BFGS-B" 法に対するグラディエント(一階(偏)微分)関数。もし NULL を指定(既定)すると、グラディエント関数は微分商を用いた数値微分で計算される。もし、比較的簡単に定義できるなら明示的に与えた方が良い。"SANN" 法に対しては、これは新しい候補パラメータを計算する関数と解釈され、もし NULL なら既定の Gaussian Markov kernel が使用される
使用する最適化手法を指示する。既定手法は "Nelder-Mead" 法
"L-BFGS-B" 法に対する変数の下限、上限を与える。既定値は lower = -Inf, upper = Inf
制御パラメータのリスト。
論理値。数値的に計算したヘッシアン(二階偏微分関数行列の値)を返すべきかどうかを指示。
fn と gr に引き渡される追加の引数(群)
Nelder-Mead 法。関数値だけを用い、頑健(例えば初期値の選択に敏感でない)であるが、相対的に遅い。微分できない関数に対してもそれなりに使える
準ニュートン法。variable metric 法とも呼ばれる。関数値とグラディエント関数を最適化関数の曲面を近似するのに使う。Broyden, Fletcher, Goldfarb, Shanno の4人が同時期に提案。
Fletcher and Reeves による共役勾配法。Polak-Ribiere と Beale-Sorenson による更新法も選択できる共役勾配法は準ニュートン法に比べると破綻しやすいが、メモリー使用量が少ないため、より大規模の最適化に使える可能性
各変数が上限・下限による制約条件を許す準ニュートン法の変種。初期値はこの制約条件を満たす必要
確率的手法であるいわゆるシミュレーテッド・アニーリング法。関数値だけを用い遅い。微分できない関数にも使える。メトロポリス法を用いる。組合せ的最適化問題(巨大な有限集合内から最適値を探す)にも使える。一般的手法ではないが、非常にラフな関数に対しても良い値を得ることができる可能性がある
非負整数。もし正なら途中結果が表示される。値が大きいほどより詳しい情報
最適化の途中で、関数 fn と gr に適用される比例定数。最適化は実際は変数 fn(par)/fnscale に対して行われる。もし負ならば、最適化は実は最大化となる
パラメータに対する比例定数(ベクトル)。最適化の際のパラメータは実際は par/parscale とされる(パラメータの変化を誇張・抑制できる)。
optim() は既定では最小化を行なう。最大化をするには制御リスト control の fnscale に負の値を与える(例えば fnscale=-1)。本来の目的関数の代わりに、値にマイナスを掛けた目的関数を用意するのもありです。
非線形連立方程式 F1(x)=F2(x)=...=Fn(x)=0 を解くには、例えば目的関数を F(x)=F1(x)^2+F2(x)^+...+Fn(x)^2 と置いて F(x) の最小化を行なう。
http://www.hppi.troitsk.ru/Kondrin/ にあるRMinpack(CRANにはまだないパッケージ)のFSolve関数を用いる。
連立方程式を解くを参照
データ d[i] にパラメータベクトル t を持つ非線形関数 f[i]=F(i, t) を当てはめるには、目的関数を SS(t)=sum((d-f)^2) と置いて t について最小化すれば良い。
目的関数をある手法で最適化した際の、返り値の適当な関数を目的関数とし、関数値のみを使用する手法である "Nelder-Mead" 法 で再帰的に最適化を行なう。例えば、非常に平坦な関数で、関数値だけではなかなか最適値が決まらないとする。最適値ではグラディエントベクトルは理論的にゼロベクトルになるべきであるから、その自乗和を目的関数に取ることが考えられる。
人生には最適化問題よりもっと大事なことがあると達観し、さっさとあきらめた方が良いかも。すべての最適化問題が解けるとは限りません。最尤推定が駄目なら、もっと簡単なモーメント型推定量を考える、等。もちろん R の optim() よりも優れたソフトがある可能性はいつでもあります。 汎用非線形最小化関数 nlm も試してみて下さい。しばしば良い結果を得るという報告があります。
既定の 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
> 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
> 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 # 最適値も確かに前より小さい
# 日本人の代表的苗字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 # 決定係数 (結構高い)
