R: nls2 错过了解决方案

Question

我正在尝试拟合双指数函数:

t = seq(0, 30, by = 0.1)
A = 20 ; B = 10 ; alpha = 0.25 ; beta = 0.01
y = A*exp(-alpha*t) + B*exp(-beta*(t))
df = as.data.frame(cbind(t,y))
ggplot(df, aes(t, y)) + geom_line() +  scale_y_continuous(limits=c(0, 50))

这个问题不能通过像log这样的简单转换来解决，所以我想使用nls2包:

library(nls2)

fo <- y ~ Ahat*exp(-alphahat*t) + Bhat*exp(-betahat*t)
fit <- nls2(fo,
            start = list(Ahat=5, Bhat=5, alphahat=0.5,betahat=0.5),
            algorithm = "brute-force",
            trace = TRUE,
            lower = c(Ahat=0, Bhat=0, alphahat=0, betahat=0),
            upper = c(Ahat=50, Bhat=50, alphahat=10,betahat=10))
fit

结果如下:

Nonlinear regression model
  model: y ~ Ahat * exp(-alphahat * t) + Bhat * exp(-betahat * t)
   data: NULL
    Ahat     Bhat alphahat  betahat 
     5.0      5.0      0.5      0.5 
 residual sum-of-squares: 37910

Number of iterations to convergence: 4 
Achieved convergence tolerance: NA

我认为我的代码有问题，因为:

• 数据:空？
• 为什么只有 4 次迭代？
• 很难想象 nls2 没有找到比起点更好的解决方案。
• 结果与解决方案相去甚远

Answer 1

根据文档，start 参数应该是定义要搜索的网格的两行的 data.frame，或者是 data.frame 包含与参数组合相对应的更多行，以测试您是否使用暴力破解。另外，nls 会给你的贴合带来麻烦，因为它是一条完美的曲线，没有噪音。暴力方法很慢，因此这里是一个减少 nls2 搜索空间的示例。然后将暴力破解 nls2 的结果用作 nls 默认算法的起始值(或者您可以使用 nls2)，添加后数据中存在一点点噪音。

## Data
t = seq(0, 30, by = 0.1)
A = 20 ; B = 10 ; alpha = 0.25 ; beta = 0.01
y = A*exp(-alpha*t) + B*exp(-beta*(t))
df = as.data.frame(cbind(t,y))

library(nls2)
fo <- y ~ Ahat*exp(-alphahat*t) + Bhat*exp(-betahat*t)

## Define the grid to search in,
## Note: decreased the grid size
grd <- data.frame(Ahat=c(10,30),
                  Bhat=c(10, 30),
                  alphahat=c(0,2),
                  betahat=c(0,1))

## Do the brute-force
fit <- nls2(fo,
            data=df,
            start = grd,
            algorithm = "brute-force",
            control=list(maxiter=100))
coef(fit)
#       Ahat       Bhat   alphahat    betahat 
# 10.0000000 23.3333333  0.0000000  0.3333333 

## Now, run through nls:
## Fails, because there is no noise
final <- nls(fo, data=df, start=as.list(coef(fit)))

## Add a little bit of noise
df$y <- df$y+rnorm(nrow(df),0,0.001)
coef((final <- nls(fo, data=df, start=as.list(coef(fit)))))
#        Ahat        Bhat    alphahat     betahat 
# 10.00034000 19.99956016  0.01000137  0.25000966 

## Plot
plot(df, col="steelblue", pch=16)
points(df$t, predict(final), col="salmon", type="l")