使用 optim 为 nls 选择初始值

Question

我在文献中看到的一种方法是使用optim()为包nls或nlme中的非线性模型选择初始值，但是，我很困惑实际实施。

举个例子，使用来自佛罗里达州 Alachua 的 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为每万人感染人数

现在，如果我想在 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() 的结果是错误的，但我看不到我的错误。感谢您的帮助。

Answer 1

主要问题是您不是 computing/returning 来自 objective 函数的平方和。然而：我觉得你真的搞反了。将 nls() 与 SSfpl 结合使用是关于 最好的 您将要进行的优化：它具有明智的启发式选择起始值（SS代表“自启动”），它为优化器提供了梯度函数。并非不可能，通过大量的工作，您可以找到更好的启发式方法来为特定系统选择起始值，但通常从 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/N-M 从良好的起始值得到合理的拟合，但显然更差；我还没分析为什么会卡在那里

使用 optim 为 nls 选择初始值

Using optim to choose initial values for nls

r

nls