如何propery指定用于optim（）或其他优化器的渐变函数
我有一个Nelder-Mead方法可以解决的优化问题，但是我也想使用BFGS或Newton-Raphson或其他需要梯度函数的方法来求解，以获得更快的速度，并希望得到更精确的估计。 我在optim / optimx文档中编写了这样一个渐变函数（我认为），但是当我将它用于BFGS我的起始值不会移动（ optim() ），否则函数彻底不运行（ optimx() ，它返回Error: Gradient function might be wrong - check it! ）。 我很抱歉有一些代码涉及再现这一点，但这里有：
这是我想要得到参数估计的函数（这是为了平滑老年死亡率，其中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)
}
你会在那里看到exp(pars) ，因为我给了log(pars)来进行优化，以便将最终的a和b限制为正数。
示例数据（1962年日本女性，如果有人很好奇）：
.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"))
以下为Nelder-Mead方法的工作：
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的矢量，相对于参数a和b 。 我还通过利用deriv()的输出得到了一个更丑陋的版本，它返回相同的答案，并且我不发布（只是为了确认衍生物是正确的）。
如果我按如下方式将其提供给optim() ，以BFGS作为方法，则估计值不会从起始值移动：
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()尝试相同的调用（未显示）时，我收到上面提到的optimx() ，并且没有返回对象。 使用"CG"指定gr = Kannisto.gr时，我会得到相同的结果。 使用"L-BFGS-B"方法，我可以得到与估计值相同的起始值，但也报告函数和梯度函数都被调用了21次，并且出现错误消息： "ERROR: BNORMAL_TERMINATION_IN_LNSRCH"
我希望能够解决这个问题的梯度函数的方式有一些小细节，因为后面的警告和optimx行为都直言不讳地暗示该函数不正确（我认为）。 我也试过maxNR()从最大化maxLik包和观察到类似的行为（初始值不动）。 任何人都可以给我一个指针？ 多谢
[编辑] @Vincent建议我与数值近似的输出进行比较：
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
如此不同的标志，并由10倍？ 我改变渐变功能以跟随：
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);
grad( function(u) KannistoLik1( c(a=u[1], b=u[2]), .Dx, .Exp ), c(1,1) );
Kannisto.gr(c(a=1,b=1), .Dx, .Exp)
这些标志是错误的：算法在朝这个方向移动时看不到任何改进，因此不会移动。
你可以使用一些计算机代数系统（这里是Maxima）为你做计算：
display2d: false;
f(a,b,x) := a * exp(b*x) / ( 1 + a * exp(b*x) );
l(a,b,d,e,x) := - d * log(f(a,b,x)) + e * f(a,b,x);
factor(diff(l(exp(a),exp(b),d,e,x),a));
factor(diff(l(exp(a),exp(b),d,e,x),b));
我只是将结果复制并粘贴到R中：
f_gradient <- function(u, .Dx, .Exp, .x.=.5:30.5) {
a <- u[1]
b <- u[1]
x <- .x.
d <- .Dx
e <- .Exp
c(
sum( (e*exp(exp(b)*x+a)-d*exp(exp(b)*x+a)-d)/(exp(exp(b)*x+a)+1)^2 ),
sum( exp(b)*x*(e*exp(exp(b)*x+a)-d*exp(exp(b)*x+a)-d)/(exp(exp(b)*x+a)+1)^2 )
)
}
library(numDeriv)
grad( function(u) KannistoLik1( c(a=u[1], b=u[2]), .Dx, .Exp ), c(1,1) )
f_gradient(c(a=1,b=1), .Dx, .Exp) # Identical
如果盲目地将梯度放在优化中，会出现数值不稳定问题：给出的解决方案是(Inf,Inf) ...为了防止它，您可以重新调整梯度（更好的解决方法是使用更少的爆炸性转换比指数，以确保参数保持正面）。
BFGSab <- optim(
log(c(a = .1, b = .1)),
fn = KannistoLik1,
gr = function(...) f_gradient(...) * 1e-3,
method = "BFGS",
.Dx = .Dx, .Exp = .Exp
)
exp(BFGSab$par) # Less precise than Nelder-Mead
上一篇: how to propery specify a gradient function for use in optim() or other optimizer