我想比较以下数据的单指数和双指数拟合。单指数函数拟合得很好,但我找不到双指数函数的起始值(我也尝试过根据合理的值手动输入一些起始值)。
有人能解释为什么双指数函数不适合这些数据以及原因吗?是否有其他可能起作用的软件包,或者这是数据的限制?
谢谢
# Dummy data
dat <- structure(list(x = c(5, 10, 30, 60, 120, 180, 240), y = c(0.000215875066736858,
0.000213143286462547, 0.000184642176655826, 0.000154069098478106,
0.000100308278926715, 6.56265655274e-05, 4.80349446283348e-05
)), class = c("tbl_df", "tbl", "data.frame"), row.names = c(NA,
-7L))
# Mono-exponential model
nls_mod_mono <- nls(y ~ statforbiology::NLS.expoDecay(x, C0, lambda1), dat)
# Predict and plot
xNew <- seq(0, 240, length.out = 100) # new grid of times
yNew <- predict(nls_mod_mono, list(x = xNew)) # predict activity
dfNew <- data.frame(x = xNew, y = yNew)
ggplot(data = dfNew, aes(x = x, y = y*100)) +
geom_line() +
geom_point(data = dat, aes(x = x, y = y*100), size = 3)
# Bi-exponential model
params <- getInitial(y ~ SSbiexp(x, intercept_initial, log_slope_initial, intercept_terminal, log_slope_terminal), data = dat)
#...failing to find initial values
nls_mod_bi <- nls(y ~ C0*exp(-lambda1*x) + C1*exp(-lambda2*x), start = list(C0 = params[1], lambda1 = exp(params[2]), C1 = params[3], lambda2 = exp(params[4])), data = dat)
我对此进行了一番研究,但没有找到它为什么不起作用的明确答案,原因显而易见:“为什么它应该很难”:
您可以通过初始化机制进行调试,或者通过添加到您的调用
debug(attr(SSbiexp, "initial"))
来跟踪其进度...内部机制使用“plinear”(部分线性)算法,它应该更加健壮,但仍然失败......trace = TRUE
getInitial
最终,我通过使用软件包中的 Nelder-Mead 优化进行强力计算,成功获得了不错的拟合效果
bbmle
...不幸的是,如果您有很多这样的数据集,并且想要以无人监督的方式对它们进行拟合,而又不需要对设置和初始参数进行大量调整,那么这将是一个相当令人头痛的问题...双指数拟合仅将对数似然提高了
logLik(fit) - logLik(fit0)
= 0.18 个单位,这个数字相当小......