如何使R语言FME包MCMC算法在非线性DDE模型拟合中收敛

问题描述 #

我使用deSolve包求解以下非线性延迟微分方程（DDE）模型：

model <- function(parms, times = seq(0, 60, 1)) {
}

已知参数包含G矩阵、U矩阵及forcing_func函数，基于观测数据与模型输出的视觉对比设定了初始参数值：

parms <- list(
  G = G,
  U = U,
  alpha_9 = 0.045,
  alpha_11 = 1.8,
  gamma = 1/3,
  tau_1 = 8,
  tau_2 = 2,
  z_e = 0.75,
  c = 2,
  i = 3,
  k_1 = 0.5,
  k_2 = 0.6,
  n = 0.1,
  q = 0.02,
  delta = 0.2)

同时拥有对应的观测数据。我尝试使用FME包的自适应延迟拒绝MCMC算法拟合部分参数，但无法实现收敛，当前代码如下：

Objective <- function(x, parset = names(x)) {
  parms[parset] <- x
  out <- model(parms) #,tout
  ## Model cost
  return(modCost(obs = data, model = out, y = "value")) 
}

Fit <- modMCMC(f = Objective,
               p = c(alpha_9 = 0.045,
                     alpha_11 = 1.8, 
                     z_e = 0.75,   
                     c = 2,       
                     i = 3,     
                     k_1 = 0.5,     
                     k_2 = 0.6,      
                     n = 0.1,     
                     q = 0.02,   
                     delta = 0.2),
               lower = c(0.01,    
                         0.01,     
                         0.01,     
                         0.1,    
                         0.1,    
                         0.1,     
                         0.1,      
                         0.01,    
                         0.001,  
                         0),     
               upper = c(1,       
                         10,     
                         1,       
                         10,      
                         10,    
                         1,       
                         1,      
                         3,     
                         10,   
                         10),  
               niter = 2000,
               updatecov = 100,
               ntrydr = 2)

收敛优化方案 #

1. 先执行参数预优化 #

直接用视觉估计的初始值跑MCMC容易陷入局部最优或收敛缓慢，建议先用modFit找到参数的最优值，以此作为MCMC的初始起点：

# 先执行参数优化，找到最优参数
opt_result <- modFit(f = Objective, 
                     p = c(alpha_9 = 0.045, alpha_11 = 1.8, z_e = 0.75, c = 2, i = 3, k_1 = 0.5, k_2 = 0.6, n = 0.1, q = 0.02, delta = 0.2),
                     lower = c(0.01, 0.01, 0.01, 0.1, 0.1, 0.1, 0.1, 0.01, 0.001, 0),
                     upper = c(1, 10, 1, 10, 10, 1, 1, 3, 10, 10))

# 使用优化后的参数作为MCMC初始值
Fit <- modMCMC(f = Objective,
               p = coef(opt_result),  # 替换为预优化的参数
               lower = c(0.01, 0.01, 0.01, 0.1, 0.1, 0.1, 0.1, 0.01, 0.001, 0),
               upper = c(1, 10, 1, 10, 10, 1, 1, 3, 10, 10),
               niter = 5000,  # 增加迭代次数，给算法足够收敛时间
               updatecov = 50,  # 更频繁更新协方差矩阵，提升自适应效率
               ntrydr = 3,  # 增加延迟拒绝尝试次数，提高样本接受率
               burnin = 1000,  # 设置燃烧期，丢弃初始未收敛的样本
               cov0 = vcov(opt_result))  # 用预优化的参数协方差作为初始 proposal 分布

2. 验证模型输出与观测数据的匹配性 #

确保model(parms)返回的结果和观测数据data的结构完全一致：

• 时间点要对应：模型输出的时间序列需覆盖观测数据的所有时间点
• 变量列名匹配：modCost的y参数指定的列名需同时存在于模型输出和观测数据中
• 检查模型是否正确实现DDE的延迟项，避免因模型逻辑错误导致输出异常，进而影响MCMC收敛

3. 优化参数拟合策略 #

• 减少参数数量：如果部分参数相关性极高（可通过pairs(opt_result)查看参数相关性），考虑固定其中部分参数，降低拟合维度
• 参数标准化：针对量级差异大的参数（如q=0.02和alpha_11=1.8），可以先做标准化处理，让参数取值范围处于同一量级，提升MCMC的采样效率

4. 监控收敛状态 #

运行MCMC后，通过以下方式验证收敛：

# 绘制轨迹图、密度图等，查看参数是否稳定
plot(Fit)

# 执行Gelman-Rubin诊断，若统计量接近1则说明收敛
gelman.diag(Fit)

若未收敛，继续增加迭代次数（如调整niter=10000），或进一步调整updatecov、ntrydr等参数。

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