警告：lmdif：info = 0。nls.lm()函数的输入参数不正确

Question

我有以下解决方案。

$$ rac{{dx}}{{dt}} = k_+(2d11 + d12 - 2x - z) - k_{-}x$$

$$ rac{{dy}}{{dt}} = k_+(2d22 + d12 - 2y - z) - k_{-}y$$

$$ rac{{dz}}{{dt}} = 2k_+(2d11 + d12 - 2x - z)(2d22 + d12 - 2y - z) - k_{-}z$$

我想用 R 编写一个程序来求解模型方程，然后将文章中表格中的数据拟合到方程的解中。我想根据 Levenberg-Marquardt 最小二乘误差算法找到最适合的 k+ 和 k−。

以下代码是我尝试这样做的。

library(deSolve)
library(minpack.lm)

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)))
  })
}

# Defining an objective function
objective = function(par){
  k_plus = par[1]
  k_minus = par[2]
  x0 = 1.71
  y0 = 1.75
  z0 = 0.62
  d11 = 1.71
  d12 = 1.75
  d22 = 0.62
  
# Vector of times
times_vec = c(1/30, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 14, 16, 18, 20)

initial_state = c(x = x0, y = y0, z = z0)
  
  # Solving the ODE
  out = lsode(y = initial_state, 
              times = times_vec,
              func = model,
              parms = c(k_plus = k_plus, k_minus = k_minus, d11 = d11, d12 = d12, d22 = d22),
              maxsteps = 100000)
  
  error = sum((out[,"x"] - xdata)^2 + (out[,"y"] - ydata)^2 + (out[,"z"] - zdata)^2)
  return(error)
}



# The table from the article
xdata = 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)
ydata = 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)
zdata = 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)


# 
fit = nls.lm(par = c(k_plus = 200,
                     k_minus = 2),
             fn = objective)

当我运行代码时，除了使用 nls.lm 函数之外，一切正常。我收到以下错误

警告：lmdif：info = 0。输入参数不正确。

我尝试过不同的参数。我查找了该函数以确保我的参数具有正确的格式，但我似乎找不到问题所在。

Answer 1

为了调试问题，我建议在尝试拟合模型之前比较模型和数据。我们可以清楚地看到，提供数值的模型结果与数据相差甚远。

library(deSolve)
library(minpack.lm)

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)))
  })
}

times_vec = c(1/30, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 14, 16, 18, 20)
k_plus = 200; k_minus = 2
x0 = 1.71; y0 = 1.75; z0 = 0.62
d11 = 1.71; d12 = 1.75; d22 = 0.62
initial_state = c(x = x0, y = y0, z = z0)

out = lsode(y = initial_state, 
            times = times_vec,
            func = model,
            parms = c(k_plus = k_plus, k_minus = k_minus, d11 = d11, d12 = d12, d22 = d22),
            maxsteps = 100000)


# The table from the article
data = data.frame(
  time = times_vec,
  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)
)

plot(out, obs=data, mfrow=c(1, 3))

因此，我建议检查：

• 模型方程和
• 数值，特别是初始状态显然与 d..-参数相同。

此外，我建议使用

lsoda

求解器而不是

lsode

。人们还可以考虑另一种优化包，例如FME，请参阅https://doi.org/10.18637/jss.v033.i03这个包是多种优化算法（包括

nls.lm

）的包装器，但对拟合ODE模型有更多支持。

解决方案

模型拟合得不太好，因为 k_plus 的起始参数太极端了。将其减少到（例如 2.0）后，技术上可以拟合模型，但曲线与数据相差太远。

然后，在包含所有参数后，就可以得到合理的拟合。

library(deSolve)
library(FME)

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)))
  })
}

times_vec <- c(1/30, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 14, 16, 18, 20)
initial_state <- c(x = 1.71, y = 1.75, z = 0.62)

# note different start value for k_plus
pstart <- c(k_plus = 1, k_minus = 1, d11 = 1.71, d12 = 1.75, d22 = 0.62)

# The table from the article
data <- data.frame(
  time = times_vec,
  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)
)

## test without fit
out <- ode(y = initial_state, times = times_vec, func = model, parms = pstart)
plot(out, obs=data, mfrow=c(1, 3))

objective <- function(par) {
  # cat(par, ", ") # uncomment for debugging
  out <- ode(y = initial_state,  times = times_vec, func = model, parms = par)
  cost <- modCost(out, data)
  # cat(cost$model, "\n") # uncomment for debugging
 cost
}

# method = "Marq" is essentially from nls.lm
# the model produces some warnings, but finally converges
# fit <- modFit(f = objective, p = pstart, method="Marq") 

# method = "Nelder" runs more smoothly in this case
fit <- modFit(f = objective, p = pstart, method="Nelder") 

# another way would be to use box constraints and then, 
# for example, L-BFGS-B or Port
#fit <- modFit(f = objective, p = pstart, 
#              lower = c(k_plus = .1, k_minus = .1, d11 = .1, d12 = .1, d22 = .1), 
#              upper = c(k_plus = 10, k_minus = 10, d11 = 10, d12 = 10, d22 = 10),
#              method="Port") 

# show detailed results
summary(fit)

# simulate model with fitted parameters
out2 <- ode(y = initial_state, times = times_vec, func = model, parms = fit$par)
plot(out2, obs=data, mfrow=c(1,3))