如何用双Logistic函数拟合跨年度时间序列？结果不稳定求解

问题背景 #

现有2016年9月1日至2020年7月30日、间隔15天的作物含水量时间序列（共96个观测值），尝试用greenbrown包的FitDoubleLogBeck函数拟合时，出现两个问题：

1. 函数无法适配整个跨年度时间序列
2. 每次运行拟合结果存在差异

问题1：无法拟合整个时间序列的解决方法 #

FitDoubleLogBeck的设计初衷是针对单年度物候周期的时序数据，默认参数逻辑仅适配一年的生长节律，跨多年的序列会让函数无法识别合理的双Logistic周期参数。

方案1：自定义双Logistic模型（推荐） #

放弃依赖FitDoubleLogBeck的年度适配逻辑，改用nls直接拟合自定义模型，用累计天数作为时间变量，适配跨年度时长：

# 生成对应日期向量
dates <- seq(as.Date("2016-09-01"), as.Date("2020-07-30"), by = "15 days")
# 转换为从起始日开始的累计天数（连续时间变量）
t <- as.numeric(dates - dates[1])

# 定义双Logistic函数（适配多年度生长-衰退周期）
double_logistic <- function(t, base_val, peak_val, rise_t, rise_rate, fall_mag, fall_t) {
  # base_val: 基线含水量；peak_val: 峰值；rise_t: 上升拐点；rise_rate: 上升速率
  # fall_mag: 下降幅度；fall_t: 下降拐点
  base_val + (peak_val - base_val)/(1 + exp(-(t - rise_t)/rise_rate)) - 
    fall_mag/(1 + exp(-(t - fall_t)/rise_rate))
}

# 基于数据范围设置合理初始值（需根据你的实际数据调整）
init_params <- list(
  base_val = min(water),
  peak_val = max(water),
  rise_t = 180,    # 第一个生长季拐点（约半年后）
  rise_rate = 30,  # 上升速率
  fall_mag = (max(water) - min(water))/2,
  fall_t = 1000    # 后期衰退拐点（约3年后）
)

# 拟合模型，放宽迭代次数限制
fit <- nls(water ~ double_logistic(t, base_val, peak_val, rise_t, rise_rate, fall_mag, fall_t),
           start = init_params,
           control = nls.control(maxiter = 1000, tol = 1e-6))

# 生成预测值并可视化
predicted <- predict(fit, newdata = data.frame(t = t))
plot(dates, water, pch = 16, col = "steelblue", main = "双Logistic拟合结果")
lines(dates, predicted, col = "firebrick", lwd = 2)

方案2：调整FitDoubleLogBeck的时间序列定义 #

如果坚持使用该函数，需将时间序列的frequency调整为总观测数对应的周期（而非年度），并手动指定初始参数：

# 用连续时间索引替代年度ts对象
x <- ts(water, start = 1, end = 96, frequency = 96)
# 手动设置初始参数（基于数据范围）
init_pars <- c(min(water), max(water), 20, 5, (max(water)-min(water))/2, 80)
# 强制使用固定初始参数，关闭随机初始化
smooth.water <- FitDoubleLogBeck(x, weighting = TRUE, hessian = FALSE, 
                                 plot = TRUE, ninit = 1, start = init_pars)

问题2：每次运行结果差异的解决方法 #

差异源于FitDoubleLogBeck的随机初始参数搜索，跨年度数据的参数空间更复杂，容易收敛到不同局部最优解。

• 增大随机初始化次数：将ninit设置为更大的值（如50或100），让函数有更多机会找到全局最优解：

  smooth.water <- FitDoubleLogBeck(x, weighting = TRUE, hessian = FALSE, 
                                   plot = TRUE, ninit = 50)
  

• 固定初始参数：如上述方案1/2中手动指定初始值，避免随机搜索的不确定性，确保每次拟合的起始点一致。
• 使用全局优化方法：如果自定义拟合，可结合nlsLM（来自minpack.lm包），它支持全局优化，减少局部最优的概率：

  library(minpack.lm)
  fit <- nlsLM(water ~ double_logistic(t, base_val, peak_val, rise_t, rise_rate, fall_mag, fall_t),
               start = init_params,
               control = nls.lm.control(maxiter = 1000))
  

内容的提问来源于stack exchange，提问作者kl40