使用nls.lm拟合ODE模型时遇lmdif参数输入错误求助
ODE参数拟合报错(lmdif: info = 0)解决方案
核心问题定位
你遇到的lmdif: info = 0. Improper input parameters错误,本质是**nls.lm要求目标函数返回残差向量,而非平方和的总和**。原代码的objective函数返回的是所有误差的平方和,不符合该函数的输入要求,这是触发报错的主要原因。
另外原代码依赖全局变量传递观测数据,虽然不报错但不够规范,也可能导致调试困难;初始参数选择如果太极端,也可能引发ODE求解的数值不稳定问题。
修正后的完整代码
library(deSolve) library(minpack.lm) # ODE模型定义 model <- function(t, state, parms) { with(as.list(c(state, parms)), { dxdt <- k_plus*(2*d11 + d12 - 2*x - z) - k_minus*x dydt <- k_plus*(2*d22 + d12 - 2*y - z) - k_minus*y dzdt <- 2*k_plus*(2*d11 + d12 - 2*x - z)*(2*d22 + d12 - 2*y - z) - k_minus*z return(list(c(dxdt, dydt, dzdt))) }) } # 修正后的目标函数:返回残差向量而非平方和 objective <- function(par, obs_data, times_vec, init_state, fixed_parms) { parms <- c(par, fixed_parms) # 求解ODE out <- lsode(y = init_state, times = times_vec, func = model, parms = parms, maxsteps = 100000) # 提取模拟值 sim_x <- out[,"x"] sim_y <- out[,"y"] sim_z <- out[,"z"] # 构造残差向量:按x、y、z的顺序拼接每个数据点的误差 residuals <- c(sim_x - obs_data$x, sim_y - obs_data$y, sim_z - obs_data$z) return(residuals) } # 观测数据整理成数据框 obs_data <- data.frame( x = c(1.71, 1.56, 1.40, 1.28, 1.14, 1.17, 1.14, 1.06, 1.10, 1.10, 1.07, 1.09, 1.06, 1.03, 1.06, 1.06, 1.05), y = c(1.75, 1.59, 1.43, 1.31, 1.17, 1.20, 1.17, 1.09, 1.13, 1.13, 1.10, 1.13, 1.09, 1.06, 1.09, 1.09, 1.08), z = c(0.62, 0.93, 1.26, 1.50, 1.78, 1.71, 1.78, 1.93, 1.86, 1.86, 1.91, 1.87, 1.93, 1.99, 1.94, 1.94, 1.95) ) # 固定参数、初始状态、时间点 fixed_parms <- c(d11 = 1.71, d12 = 1.75, d22 = 0.62) init_state <- c(x = 1.71, y = 1.75, z = 0.62) times_vec <- c(1/30, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 14, 16, 18, 20) # 调用nls.lm拟合:调整初始参数为更合理的值 fit <- nls.lm(par = c(k_plus = 10, k_minus = 0.5), fn = objective, obs_data = obs_data, times_vec = times_vec, init_state = init_state, fixed_parms = fixed_parms) # 查看拟合结果 summary(fit)
关键修改说明
- 目标函数返回值调整:把原来的平方和总和改成直接返回每个观测点与模拟点的残差组成的向量,这是
nls.lm的强制要求——函数需要返回长度等于所有观测数据点总数的残差向量(这里是17*3=51个残差)。 - 参数传递规范化:将观测数据、时间点、初始状态、固定参数都作为参数传入目标函数,不再依赖全局变量,代码更健壮。
- 初始参数优化:把k_plus的初始值从200下调到10,避免ODE求解时出现数值爆炸导致的不稳定问题,你可以根据拟合结果进一步调整初始值。
- 数据结构优化:把观测数据整理成数据框,更便于管理和传递。
运行修正后的代码后,应该能正常执行拟合,你可以通过summary(fit)查看最优的k_plus和k_minus参数值,以及拟合的统计信息。
内容的提问来源于stack exchange,提问作者Nick
相关产品推荐
相关产品推荐

