使用optim为nls选择初始值

我在文献中所看到的一种方法是使用optim()选择在包非线性模型的初始值nlsnlme,但是,我通过实际执行不解。

以佛罗里达州阿拉楚阿的 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切换nlsoptim+ 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 从好的起始值得到合理的拟合,但显然更糟;我还没有分析它为什么卡在那里。

以上是使用optim为nls选择初始值的全部内容。
THE END
分享
二维码
< <上一篇
下一篇>>