使用optim为nls选择初始值
我在文献中所看到的一种方法是使用optim()选择在包非线性模型的初始值nls或nlme,但是,我通过实际执行不解。
以佛罗里达州阿拉楚阿的 COVID 数据为例:
dat=data.frame(x=seq(1,10,1), y=c(27.9,23.1,24.6,33.0,48.0,136.4,243.4,396.7,519.9,602.8))
x 是时间点,y 是每 10,000 人感染的人数
现在,如果我想在 nls 中拟合一个四参数逻辑模型,我可以使用
n1 <- nls(y ~ SSfpl(x, A, B, M, S), data = dat)
但是现在想象一下参数估计对初始值高度敏感,所以我想优化我的方法。这将如何实现?
我想尝试的方法如下
fun_to_optim <- function(data, guess){
x = data$x
y = data$y
A = guess[1]
B = guess[2]
M = guess[3]
S = guess[4]
y = A + (B-A)/(1+exp((M-x)/S))
return(-sum(y)) }
optim(fn=fun_to_optim, data=dat,
par=c(10,10,10,10),
method="Nelder-Mead")
结果optim()是错误的,但我看不到我的错误。感谢您提供任何帮助。
回答
主要问题是您没有从目标函数计算/返回平方和。 但是:我认为你真的倒退了。使用nls()withSSfpl是您在优化方面要做的最好的事情:它具有用于选择起始值的明智启发式(SS代表“自启动”),并且它为优化器提供了梯度函数。并非不可能,通过大量工作,您可以找到更好的启发式方法来为特定系统选择起始值,但通常从+ Nelder-Mead切换nls到optim+ Nelder-Mead 会使您的情况比开始时更糟(下图)。
fun_to_optim <- function(data, guess){
x = data$x
y = data$y
A = guess[1]
B = guess[2]
M = guess[3]
S = guess[4]
y_pred = A + (B-A)/(1+exp((M-x)/S))
return(sum((y-y_pred)^2))
}
符合optim()第(1)你的建议的起始值; (2) 更接近正确值的更好的起始值(你可以通过了解函数的几何来获得这些值中的大部分——例如A是左渐近线,B是右渐近线,M是中点,S是尺度);(3) 与 #2 相同,但使用 BFGS 而不是 Nelder-Mead。
opt1 <- optim(fn=fun_to_optim, data=dat,
par=c(A=10,B=10,M=10,S=10),
method="Nelder-Mead")
opt2 <- optim(fn=fun_to_optim, data=dat,
par=c(A=10,B=500,M=10,S=1),
method = "Nelder-Mead")
opt3 <- optim(fn=fun_to_optim, data=dat,
par=c(A=10,B=500,M=10,S=1),
method = "BFGS")
结果:
xvec <- seq(1,10,length=101)
plot(y~x, data=dat)
lines(xvec, predict(n1, newdata=data.frame(x=xvec)))
p1 <- with(as.list(opt1$par), A + (B-A)/(1+exp((M-xvec)/S)))
lines(xvec, p1, col=2)
p2 <- with(as.list(opt2$par), A + (B-A)/(1+exp((M-xvec)/S)))
lines(xvec, p2, col=4)
p3 <- with(as.list(opt3$par), A + (B-A)/(1+exp((M-xvec)/S)))
lines(xvec, p3, col=6)
legend("topleft", col=c(1,2,4,6), lty=1,
legend=c("nls","NM (bad start)", "NM", "BFGS"))
nls和良好的起始值 + BFGS 重叠,并提供良好的拟合optim/Nelder-Mead 来自糟糕的起始值是绝对可怕的——收敛在一条恒定线上optim/NM 从好的起始值得到合理的拟合,但显然更糟;我还没有分析它为什么卡在那里。