関数の最大・最小化
をテンプレートにして作成
[
トップ
] [
新規
|
一覧
|
検索
|
最終更新
|
ヘルプ
]
開始行:
//間瀬茂
COLOR(red){SIZE(30){汎用最適化関数 optim() の使用法}}
-統計学では、最尤推定や最小自乗法等、ある目的関数の最適化...
しばしば起こる。COLOR(red){optim()} は R の汎用最適化関数...
-一変数関数専用の最適化関数 COLOR(red){optimize} もある。
-線形不等式制約の場合は [[線形不等式制約付きの最適化関数 ...
-同種の関数として [[汎用非線形最小化関数 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"...
lower = -Inf, upper = Inf,
control = list(), hessian = FALSE, ...)
COLOR(red){optim()} は再帰的に使うことができる(目的関数自...
**引数の意味 [#w4e72c3a]
*** par [#yc3cb4c4]
目的関数のベクトル引数に対する初期値。この選択は一般に試...
*** fn [#o8af719e]
最小化(既定)もしくは最大化すべき目的変数。最初の(ベクトル...
*** gr [#bc97c2ae]
"BFGS", "CG", "L-BFGS-B" 法に対するグラディエント(一階(偏...
*** method [#sf1e857d]
使用する最適化手法を指示する。既定手法は COLOR(magenta){"...
*** lower, upper [#h6b41a03]
COLOR(magenta){"L-BFGS-B" 法}に対する変数の下限、上限を与...
*** 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 m...
*** "CG" 法 [#k48c0acb]
Fletcher and Reeves によるCOLOR(magenta){共役勾配法}。Pol...
*** "L-BFGS-B" 法 [#c303cc6e]
各変数が上限・下限によるCOLOR(magenta){制約条件を許す準ニ...
*** "SANN" 法 [#z688fed4]
確率的手法であるいわゆるCOLOR(magenta){シミュレーテッド・...
** 制御パラメータリスト control の意味 [#uaacef6e]
*** trace [#j62e2900]
非負整数。もし正なら途中結果が表示される。値が大きいほど...
*** fnscal [#gf14ec83]
最適化の途中で、関数 COLOR(red){fn} と COLOR(red){gr} に...
*** parscale [#h3018f50]
パラメータに対する比例定数(ベクトル)。最適化の際のパラメ...
* 注意 [#e28ecc95]
** 最大化するには？ [#ecf7726d]
SIZE(20){optim()} は既定では最小化を行なう。最大化をする...
** 非線形連立方程式を解くには？ [#h9b8656f]
*** optim関数 [#hb12e05f]
非線形連立方程式 COLOR(red){F1(x)=F2(x)=...=Fn(x)=0} を解...
*** RMinpackパッケージ [#s8a530f2]
http://www.hppi.troitsk.ru/Kondrin/ にあるRMinpack(CRANに...
*** nleqslvパッケージ [#v7f2c799]
[[連立方程式を解く]]を参照
** 非線形回帰をするには? [#gde33b2f]
データ COLOR(red){d[i]} にパラメータベクトル COLOR(red){t...
COLOR(red){SS(t)=sum((d-f)^2)} と置いて COLOR(red){t} に...
** 変数に制約がある時は? [#i1f12672]
-変数に制約がある時は最適化は一層困難になります。矩形型の...
-しばしば使われる方法は、制約条件を表す関数 COLOR(red){g(...
-線形不等式制約の場合は [[線形不等式制約付きの最適化関数 ...
**本当に最適値が求まったの? [#ybf88577]
-これは意外に難しい点です。如何なる数値的最適化手法も、CO...
-さらに計算に於けるCOLOR(magenta){数値誤差の問題}も無視で...
-普通COLOR(magneta){最適値があると信じ}て計算するのでしょ...
-また、本当に真の最適値が必要なのでしょうか(良くある例で...
-まず収束したかどうかを結果の収束判定コード COLOR(red){co...
-収束判定用の基準値 COLOR(red){abstol}, COLOR(red){reltol...
-もし、収束しないままに打ち切られていたら、たとえ、尤もら...
-微分可能な関数ならば、(制約付き最適化は例外かもしれない)...
-いろいろ書くとやる気を無くしそうですが、心配しないで下さ...
-統計学ではランダムなデータを扱いますから、誤差が大きけれ...
** 初期値パラメータの選び方は？ [#kbc1995a]
***神頼み、人頼み [#d73dac7a]
-基本的に試行錯誤を繰り返して選ぶより仕方がないと心得る
-問題の性質からありそうな範囲の大雑把な見当をつけるべきで...
-教科書には書いてないお勧めの方法として、COLOR(magenta){...
***システマティックな探索 [#ae232611]
-パラメータの範囲を適当に等分して、各升目毎に、目的関数の...
-1変数ならば COLOR(red){plot()} 関数でグラフを描いてみる。
-2変数ならば、計算結果を行列に表現した上で、 COLOR(red){i...
*** 多段階最適化 [#n455d54b]
-COLOR(magenta){"Nelder-Mead" 法} や COLOR(magenta){"SANN...
-次に、それを初期値として、改めて他の切れ味の鋭い(ただし...
-微分できないような関数ではCOLOR(magenta){"Nelder-Mead" ...
*** それでも駄目なら (1) [#ta9ecccd]
-問題の設定を変えてみる。問題を同値なより単純な問題に整理...
-できれば、変数の数を減らした問題に置き換える。
-変数を置き換えて制約条件を除く。例えば、確率パラメータ C...
-少しは頭を使い、実は紙と鉛筆で解ける問題ではないか反省し...
-繰り返し回数を増やしてみる。
-COLOR(red){parscale} を大きくし、一回毎のパラメータの変...
-COLOR(red){fnscale} を小さな(正)数にし、関数値の凹凸を誇...
*** それでも駄目なら (2) [#h8e5e6fb]
目的関数をある手法で最適化した際の、COLOR(red){返り値の適...
-optim(par, fn=optim(par, fn=SS)$value) # ただし非常に時...
*** それでも駄目なら (3) [#s265c525]
人生には最適化問題よりもっと大事なことがあると達観し、さ...
[[汎用非線形最小化関数 nlm]] も試してみて下さい。しばしば...
*実例 [#df5368fe]
** 実例 (1) 負の二項分布のパラメータの最尤推定 [#s8619974]
既定の Nelder-Mead 法で対数尤度を最大化(オプション fnscal...
d <- c(1612,164,71,47,28,17,12,12,5,7,6,3,3,13) # データ...
NB <- function (s, p) { # パラメータ (size, prob)=(s, p)...
P <- dnbinom(0:12, size=s, prob=p, log=T) # 個数 0,1,....
P[14] <- pnbinom(12, size=s, prob=p, lower.tail=F, log...
return(P)
}
LL <- function (x) return(sum(NB(x[1],x[2])*d)) # 目的関...
optim(c(0.2, 0.5), LL, control=list(fnscale=-1)) # 初期...
$par
[1] 0.1155133 0.1553361 # (size, prob) パラメータの最尤...
$value
[1] -1705.324 # その時の対数尤度値
$counts
function gradient
63 NA # Nelder-Mead 法繰り返し数63回...
$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 2...
[7] 14.546131 10.734116 8.064296 6.142199 4...
[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 * (x...
}
# 以下各種手法による最適パラメータと最適値を紹介
# 多数決により、最適値は 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") # "BFG...
$par
[1] 1 1
$value
[1] 9.594955e-18
> optim(c(-1.2, 1), fr, NULL, method = "BFGS", hessian =...
$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 = li...
$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(m...
temp = 20, parscale = 20)) # "SA...
> 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位の頻度...
# 日本人の代表的苗字100位(電子電話帳を集計したデータ。単...
d <- c(456430,403506,335288,314770,256706,255876,254662,...
197460,193503,169617,152065,149006,143552,137475,1...
114802,110430,108369,108345,105778,102647, 97704, ...
90925, 89856, 89818, 87815, 86992, 86695, 86234, ...
75826, 75264, 74510, 74352, 73185, 72640, 72569, ...
70082, 69904, 68661, 67852, 67571, 65830, 64234, ...
59967, 58141, 57568, 57037, 56651, 56538, 56324, ...
53284, 52858, 50891, 50499, 50349, 50180, 49474, ...
48724, 48386, 48329, 47744, 47094, 46923, 46858, ...
45682, 45164, 44731, 44650, 44641, 44222, 43908, ...
d <- d/sum(d) # 頻度
zipf <- function (x) { # Zipf 分布 p(i) = c/(i^x), i=1,2...
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, parscal...
> res$par
[1] 1.030962 # 非線形最小自乗推定値
> 1-SS(res$par)/sum(d^2)
[1] 0.9336638 # 決定係数 (結構高い)
終了行:
//間瀬茂
COLOR(red){SIZE(30){汎用最適化関数 optim() の使用法}}
-統計学では、最尤推定や最小自乗法等、ある目的関数の最適化...
しばしば起こる。COLOR(red){optim()} は R の汎用最適化関数...
-一変数関数専用の最適化関数 COLOR(red){optimize} もある。
-線形不等式制約の場合は [[線形不等式制約付きの最適化関数 ...
-同種の関数として [[汎用非線形最小化関数 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"...
lower = -Inf, upper = Inf,
control = list(), hessian = FALSE, ...)
COLOR(red){optim()} は再帰的に使うことができる(目的関数自...
**引数の意味 [#w4e72c3a]
*** par [#yc3cb4c4]
目的関数のベクトル引数に対する初期値。この選択は一般に試...
*** fn [#o8af719e]
最小化(既定)もしくは最大化すべき目的変数。最初の(ベクトル...
*** gr [#bc97c2ae]
"BFGS", "CG", "L-BFGS-B" 法に対するグラディエント(一階(偏...
*** method [#sf1e857d]
使用する最適化手法を指示する。既定手法は COLOR(magenta){"...
*** lower, upper [#h6b41a03]
COLOR(magenta){"L-BFGS-B" 法}に対する変数の下限、上限を与...
*** 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 m...
*** "CG" 法 [#k48c0acb]
Fletcher and Reeves によるCOLOR(magenta){共役勾配法}。Pol...
*** "L-BFGS-B" 法 [#c303cc6e]
各変数が上限・下限によるCOLOR(magenta){制約条件を許す準ニ...
*** "SANN" 法 [#z688fed4]
確率的手法であるいわゆるCOLOR(magenta){シミュレーテッド・...
** 制御パラメータリスト control の意味 [#uaacef6e]
*** trace [#j62e2900]
非負整数。もし正なら途中結果が表示される。値が大きいほど...
*** fnscal [#gf14ec83]
最適化の途中で、関数 COLOR(red){fn} と COLOR(red){gr} に...
*** parscale [#h3018f50]
パラメータに対する比例定数(ベクトル)。最適化の際のパラメ...
* 注意 [#e28ecc95]
** 最大化するには？ [#ecf7726d]
SIZE(20){optim()} は既定では最小化を行なう。最大化をする...
** 非線形連立方程式を解くには？ [#h9b8656f]
*** optim関数 [#hb12e05f]
非線形連立方程式 COLOR(red){F1(x)=F2(x)=...=Fn(x)=0} を解...
*** RMinpackパッケージ [#s8a530f2]
http://www.hppi.troitsk.ru/Kondrin/ にあるRMinpack(CRANに...
*** nleqslvパッケージ [#v7f2c799]
[[連立方程式を解く]]を参照
** 非線形回帰をするには? [#gde33b2f]
データ COLOR(red){d[i]} にパラメータベクトル COLOR(red){t...
COLOR(red){SS(t)=sum((d-f)^2)} と置いて COLOR(red){t} に...
** 変数に制約がある時は? [#i1f12672]
-変数に制約がある時は最適化は一層困難になります。矩形型の...
-しばしば使われる方法は、制約条件を表す関数 COLOR(red){g(...
-線形不等式制約の場合は [[線形不等式制約付きの最適化関数 ...
**本当に最適値が求まったの? [#ybf88577]
-これは意外に難しい点です。如何なる数値的最適化手法も、CO...
-さらに計算に於けるCOLOR(magenta){数値誤差の問題}も無視で...
-普通COLOR(magneta){最適値があると信じ}て計算するのでしょ...
-また、本当に真の最適値が必要なのでしょうか(良くある例で...
-まず収束したかどうかを結果の収束判定コード COLOR(red){co...
-収束判定用の基準値 COLOR(red){abstol}, COLOR(red){reltol...
-もし、収束しないままに打ち切られていたら、たとえ、尤もら...
-微分可能な関数ならば、(制約付き最適化は例外かもしれない)...
-いろいろ書くとやる気を無くしそうですが、心配しないで下さ...
-統計学ではランダムなデータを扱いますから、誤差が大きけれ...
** 初期値パラメータの選び方は？ [#kbc1995a]
***神頼み、人頼み [#d73dac7a]
-基本的に試行錯誤を繰り返して選ぶより仕方がないと心得る
-問題の性質からありそうな範囲の大雑把な見当をつけるべきで...
-教科書には書いてないお勧めの方法として、COLOR(magenta){...
***システマティックな探索 [#ae232611]
-パラメータの範囲を適当に等分して、各升目毎に、目的関数の...
-1変数ならば COLOR(red){plot()} 関数でグラフを描いてみる。
-2変数ならば、計算結果を行列に表現した上で、 COLOR(red){i...
*** 多段階最適化 [#n455d54b]
-COLOR(magenta){"Nelder-Mead" 法} や COLOR(magenta){"SANN...
-次に、それを初期値として、改めて他の切れ味の鋭い(ただし...
-微分できないような関数ではCOLOR(magenta){"Nelder-Mead" ...
*** それでも駄目なら (1) [#ta9ecccd]
-問題の設定を変えてみる。問題を同値なより単純な問題に整理...
-できれば、変数の数を減らした問題に置き換える。
-変数を置き換えて制約条件を除く。例えば、確率パラメータ C...
-少しは頭を使い、実は紙と鉛筆で解ける問題ではないか反省し...
-繰り返し回数を増やしてみる。
-COLOR(red){parscale} を大きくし、一回毎のパラメータの変...
-COLOR(red){fnscale} を小さな(正)数にし、関数値の凹凸を誇...
*** それでも駄目なら (2) [#h8e5e6fb]
目的関数をある手法で最適化した際の、COLOR(red){返り値の適...
-optim(par, fn=optim(par, fn=SS)$value) # ただし非常に時...
*** それでも駄目なら (3) [#s265c525]
人生には最適化問題よりもっと大事なことがあると達観し、さ...
[[汎用非線形最小化関数 nlm]] も試してみて下さい。しばしば...
*実例 [#df5368fe]
** 実例 (1) 負の二項分布のパラメータの最尤推定 [#s8619974]
既定の Nelder-Mead 法で対数尤度を最大化(オプション fnscal...
d <- c(1612,164,71,47,28,17,12,12,5,7,6,3,3,13) # データ...
NB <- function (s, p) { # パラメータ (size, prob)=(s, p)...
P <- dnbinom(0:12, size=s, prob=p, log=T) # 個数 0,1,....
P[14] <- pnbinom(12, size=s, prob=p, lower.tail=F, log...
return(P)
}
LL <- function (x) return(sum(NB(x[1],x[2])*d)) # 目的関...
optim(c(0.2, 0.5), LL, control=list(fnscale=-1)) # 初期...
$par
[1] 0.1155133 0.1553361 # (size, prob) パラメータの最尤...
$value
[1] -1705.324 # その時の対数尤度値
$counts
function gradient
63 NA # Nelder-Mead 法繰り返し数63回...
$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 2...
[7] 14.546131 10.734116 8.064296 6.142199 4...
[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 * (x...
}
# 以下各種手法による最適パラメータと最適値を紹介
# 多数決により、最適値は 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") # "BFG...
$par
[1] 1 1
$value
[1] 9.594955e-18
> optim(c(-1.2, 1), fr, NULL, method = "BFGS", hessian =...
$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 = li...
$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(m...
temp = 20, parscale = 20)) # "SA...
> 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位の頻度...
# 日本人の代表的苗字100位(電子電話帳を集計したデータ。単...
d <- c(456430,403506,335288,314770,256706,255876,254662,...
197460,193503,169617,152065,149006,143552,137475,1...
114802,110430,108369,108345,105778,102647, 97704, ...
90925, 89856, 89818, 87815, 86992, 86695, 86234, ...
75826, 75264, 74510, 74352, 73185, 72640, 72569, ...
70082, 69904, 68661, 67852, 67571, 65830, 64234, ...
59967, 58141, 57568, 57037, 56651, 56538, 56324, ...
53284, 52858, 50891, 50499, 50349, 50180, 49474, ...
48724, 48386, 48329, 47744, 47094, 46923, 46858, ...
45682, 45164, 44731, 44650, 44641, 44222, 43908, ...
d <- d/sum(d) # 頻度
zipf <- function (x) { # Zipf 分布 p(i) = c/(i^x), i=1,2...
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, parscal...
> res$par
[1] 1.030962 # 非線形最小自乗推定値
> 1-SS(res$par)/sum(d^2)
[1] 0.9336638 # 決定係数 (結構高い)
ページ名:
![[RjpWiki]](image/../Rlogo.png "[RjpWiki]")
