Nelder-Meadメソッドが解決する最適化の問題がありますが、BFGSまたはNewton-Raphsonなどを使用して解決したいと思います。より速いスピードのための勾配関数、そしてより正確な見積もりが望まれます。私はoptim/optimxのドキュメントの例に従って（私は思った）そのような勾配関数を書きましたが、BFGSと一緒に使用すると、私の開始値が動かず（optim()）、そうでなければ関数が実行されません。optimx() 、それはError: Gradient function might be wrong - check it!を返す）。optim（）やその他のオプティマイザで使用するグラディエント関数を指定する方法
これはパラメータ推定値を取得したい機能です（これは、老齢死亡率を平滑化するための関数です。xは年齢、 ）80歳から始まる：
KannistoMu <- function(pars, x = .5:30.5){
a <- pars["a"]
b <- pars["b"]
(a * exp(b * x))/(1 + a * exp(b * x))
}
そして、ここで観測された率（死亡、露出オーバー.Dx、.Expと定義）からそれを推定するための対数尤度関数です：
KannistoLik1 <- function(pars, .Dx, .Exp, .x. = .5:30.5){
mu <- KannistoMu(exp(pars), x = .x.)
# take negative and minimize it (default optimizer behavior)
-sum(.Dx * log(mu) - .Exp * mu, na.rm = TRUE)
}
あなたがそこBECAでexp(pars)を見ます最終的にaとbが肯定的になるように、私はlog(pars)を使って最適化を行います。
例データ（1962日本の女性、誰でも好奇心旺盛であれば）：Nelder-Mead方法について
.Dx <- structure(c(10036.12, 9629.12, 8810.11, 8556.1, 7593.1, 6975.08,
6045.08, 4980.06, 4246.06, 3334.04, 2416.03, 1676.02, 1327.02,
980.02, 709, 432, 350, 217, 134, 56, 24, 21, 10, 8, 3, 1, 2,
1, 0, 0, 0), .Names = c("80", "81", "82", "83", "84", "85", "86",
"87", "88", "89", "90", "91", "92", "93", "94", "95", "96", "97",
"98", "99", "100", "101", "102", "103", "104", "105", "106",
"107", "108", "109", "110"))
.Exp <- structure(c(85476.0333333333, 74002.0866666667, 63027.5183333333,
53756.8983333333, 44270.9, 36749.85, 29024.9333333333, 21811.07,
16912.315, 11917.9583333333, 7899.33833333333, 5417.67, 3743.67833333333,
2722.435, 1758.95, 1043.985, 705.49, 443.818333333333, 223.828333333333,
93.8233333333333, 53.1566666666667, 27.3333333333333, 16.1666666666667,
10.5, 4.33333333333333, 3.16666666666667, 3, 2.16666666666667,
1.5, 0, 1), .Names = c("80", "81", "82", "83", "84", "85", "86",
"87", "88", "89", "90", "91", "92", "93", "94", "95", "96", "97",
"98", "99", "100", "101", "102", "103", "104", "105", "106",
"107", "108", "109", "110"))
次作品：
：NMab <- optim(log(c(a = .1, b = .1)),
fn = KannistoLik1, method = "Nelder-Mead",
.Dx = .Dx, .Exp = .Exp)
exp(NMab$par)
# these are reasonable estimates
a b
0.1243144 0.1163926
これは私が思いついた勾配関数であります
Kannisto.gr <- function(pars, .Dx, .Exp, x = .5:30.5){
a <- exp(pars["a"])
b <- exp(pars["b"])
d.a <- (a * exp(b * x) * .Exp + (-a * exp(b * x) - 1) * .Dx)/
(a^3 * exp(2 * b * x) + 2 * a^2 * exp(b * x) + a)
d.b <- (a * x * exp(b * x) * .Exp + (-a * x * exp(b * x) - x) * .Dx)/
(a^2 * exp(2 * b * x) + 2 * a * exp(b * x) + 1)
-colSums(cbind(a = d.a, b = d.b), na.rm = TRUE)
}
出力は長さ2のベクトルであり、pアラムターaおよびb。私はまた、同じ回答を返すderiv()の出力を利用して醜いバージョンが届いていますが、投稿していません（派生商品が正しいことを確認するだけです）。
方法としてBFGSで、次のように私はoptim()に供給した場合、推定値は開始値から移動しない：
BFGSab <- optim(log(c(a = .1, b = .1)),
fn = KannistoLik1, gr = Kannisto.gr, method = "BFGS",
.Dx = .Dx, .Exp = .Exp)
# estimates do not change from starting values:
exp(BFGSab$par)
a b
0.1 0.1
私は、出力の$counts要素を見て、それがいることを言いますKannistoLik1()は31回、Kannisto.gr()は1回だけ呼び出されました。 $convergenceは0なので、それは収束していると思います。私は公差などを減らし、何も変わらない。同じ電話をoptimx()（表示されていません）で試してみると、私は上記のことを受け取り、何も返されません。 gr = Kannisto.grを"CG"と指定しても同じ結果が得られます。"L-BFGS-B"方法で、私は戻って見積もりと同じ開始値を取得し、それはまた、関数と勾配の両方が21回呼び出されたことが報告され、エラーメッセージがあります： "ERROR: BNORMAL_TERMINATION_IN_LNSRCH"
私はいくつか存在することを望んでいるがこの後の警告とoptimxのふるまいは、関数が単純に正しいとは思っていない（私は思う）ように、これを解決するグラディエント関数が書かれている方法の細かいディテールです。また、maxLikパッケージからmaxNR()最大化を試みましたが、同様の動作が観察されました（開始値は移動しません）。誰も私にポインタを与えることができますか？多くの義務
[編集] @Vincent私は数値近似からの出力と比較提案：、及び10倍オフ
library(numDeriv)
grad(function(u) KannistoLik1(c(a=u[1], b=u[2]), .Dx, .Exp), log(c(.1,.1)))
[1] -14477.40 -7458.34
Kannisto.gr(log(c(a=.1,b=.1)), .Dx, .Exp)
a b
144774.0 74583.4
ので、異なる符号を？私は追随する勾配関数を変更します。
Kannisto.gr2 <- function(pars, .Dx, .Exp, x = .5:30.5){
a <- exp(pars["a"])
b <- exp(pars["b"])
d.a <- (a * exp(b * x) * .Exp + (-a * exp(b * x) - 1) * .Dx)/
(a^3 * exp(2 * b * x) + 2 * a^2 * exp(b * x) + a)
d.b <- (a * x * exp(b * x) * .Exp + (-a * x * exp(b * x) - x) * .Dx)/
(a^2 * exp(2 * b * x) + 2 * a * exp(b * x) + 1)
colSums(cbind(a=d.a,b=d.b), na.rm = TRUE)/10
}
Kannisto.gr2(log(c(a=.1,b=.1)), .Dx, .Exp)
# same as numerical:
a b
-14477.40 -7458.34
オプティマイザでそれを試してみてください。
BFGSab <- optim(log(c(a = .1, b = .1)),
fn = KannistoLik1, gr = Kannisto.gr2, method = "BFGS",
.Dx = .Dx, .Exp = .Exp)
# not reasonable results:
exp(BFGSab$par)
a b
Inf Inf
# and in fact, when not exp()'d, they look oddly familiar:
BFGSab$par
a b
-14477.40 -7458.34
ヴィンセントの答えの後、私は勾配関数を再スケーリングし、正のパラメータを維持するためにabs()の代わりexp()を使用しました。最新の、そしてより良い行い目的と勾配機能は：
KannistoLik2 <- function(pars, .Dx, .Exp, .x. = .5:30.5){
mu <- KannistoMu.c(abs(pars), x = .x.)
# take negative and minimize it (default optimizer behavior)
-sum(.Dx * log(mu) - .Exp * mu, na.rm = TRUE)
}
# gradient, to be down-scaled in `optim()` call
Kannisto.gr3 <- function(pars, .Dx, .Exp, x = .5:30.5){
a <- abs(pars["a"])
b <- abs(pars["b"])
d.a <- (a * exp(b * x) * .Exp + (-a * exp(b * x) - 1) * .Dx)/
(a^3 * exp(2 * b * x) + 2 * a^2 * exp(b * x) + a)
d.b <- (a * x * exp(b * x) * .Exp + (-a * x * exp(b * x) - x) * .Dx)/
(a^2 * exp(2 * b * x) + 2 * a * exp(b * x) + 1)
colSums(cbind(a = d.a, b = d.b), na.rm = TRUE)
}
# try it out:
BFGSab2 <- optim(
c(a = .1, b = .1),
fn = KannistoLik2,
gr = function(...) Kannisto.gr3(...) * 1e-7,
method = "BFGS",
.Dx = .Dx, .Exp = .Exp
)
# reasonable:
BFGSab2$par
a b
0.1243249 0.1163924
# better:
KannistoLik2(exp(NMab1$par),.Dx = .Dx, .Exp = .Exp) > KannistoLik2(BFGSab2$par,.Dx = .Dx, .Exp = .Exp)
[1] TRUE
これは私が期待していたよりもはるかに速く解決した、と私はカップルのトリックよりも多くのことを学びました。おかげさまでヴィンセント！
グラジエントが正しいかどうかを確認するには、数値の近似値、たとえば 'library（numDeriv）;を比較できます。 （1）、b = u [2]）、.Dx、.Exp）、c（1,1））;勾配関数（関数（u）KannistoLik1 Kannisto.gr（c（a = 1、b = 1）、.Dx、.Exp） '。兆候は間違っています。アルゴリズムはこの方向に動いたときに改善が見られないため、動きません。 –
Vincentに感謝します。それを試して、上記の結果を投稿します –