我有一个优化问题,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)
}
你在那里贝科看到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的向量,相对于p的变化参数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
元素,它说,被称为31次,Kannisto.gr()
只有1次。 $convergence
是0
,所以我想它认为它收敛了(如果我给予不太合理的开始,他们也保持放置)。我减少了宽容等,没有任何变化。当我在optimx()
(未显示)尝试同样的呼叫时,我收到上面提到的战斗,并且没有任何对象返回。指定gr = Kannisto.gr
与"CG"
时,我会得到相同的结果。随着"L-BFGS-B"
方法我会得到相同的初始值追溯到估计,但也有报道,这两个函数和梯度被称为21次,并有一条错误消息: "ERROR: BNORMAL_TERMINATION_IN_LNSRCH"
我希望有一些稍微详细一点的梯度函数的写法,将解决这个问题,因为这个后来的警告和optimx
行为都直言不讳地暗示该函数根本不对(我认为)。我也尝试了maxLik
包中的maxNR()
最大化器,并观察到类似的行为(起始值不移动)。任何人都可以给我一个指针?非常感谢
[编辑] @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); (函数(u)KannistoLik1(c(a = u [1],b = u [2]),.Dx,.Exp),c(1,1)); Kannisto.gr(c(a = 1,b = 1),.Dx,.Exp)'。这些标志是错误的:算法在朝这个方向移动时看不到任何改进,因此不会移动。 – 2012-07-24 01:57:21
谢谢文森特。试了一下,会发布结果上面 – 2012-07-24 02:04:15