使用 nls.lm 拟合 ODE 模型的参数和初始条件

Question

我目前正在尝试按照此处的教程 (http://www.r-bloggers.com/learning-r-parameter-fitting-for-models-involving-differential-equations/) 使用 pkg-minpack.lm 中的 Levenberg-Marquardt 例程 (nls.lm) 拟合 ODE 函数响应。

在示例中，他首先设置了一个函数 rxnrate 来拟合数据，我对其进行了如下修改：

library(ggplot2) #library for plotting
library(reshape2) # library for reshaping data (tall-narrow <-> short-wide)
library(deSolve) # library for solving differential equations
library(minpack.lm) # library for least squares fit using levenberg-marquart algorithm
# prediction of concentration
# rate function
rxnrate=function(t,c,parms){
  # rate constant passed through a list called parms
  k1=parms$k1
  k2=parms$k2
  k3=parms$k3

 # c is the concentration of species

 # derivatives dc/dt are computed below
  r=rep(0,length(c))
  r[1]=-k1*c["A"]  #dcA/dt
  r[2]=k1*c["A"]-k2*c["B"]+k3*c["C"] #dcB/dt
  r[3]=k2*c["B"]-k3*c["C"] #dcC/dt

  # the computed derivatives are returned as a list
  # order of derivatives needs to be the same as the order of species in c
  return(list(r))

}  

我的问题是每个状态的初始条件也可以视为估计参数。但是，它目前无法正常工作。 以下是我的代码：

# function that calculates residual sum of squares
ssq=function(myparms){

  # inital concentration 
  cinit=c(A=myparms[4],B=0,C=0)

  # time points for which conc is reported
  # include the points where data is available
  t=c(seq(0,5,0.1),df$time)
  t=sort(unique(t))
  # parms from the parameter estimation routine
  k1=myparms[1]
  k2=myparms[2]
  k3=myparms[3]
  # solve ODE for a given set of parameters
   out=ode(y=cinit,times=t,func=rxnrate,parms=list(k1=k1,k2=k2,k3=k3))


  # Filter data that contains time points where data is available
  outdf=data.frame(out)
  outdf=outdf[outdf$time %in% df$time,]
  # Evaluate predicted vs experimental residual
  preddf=melt(outdf,id.var="time",variable.name="species",value.name="conc")
  expdf=melt(df,id.var="time",variable.name="species",value.name="conc")
  ssqres=preddf$conc-expdf$conc

  # return predicted vs experimental residual
  return(ssqres)

}

# parameter fitting using levenberg marquart algorithm
# initial guess for parameters
 myparms=c(k1=0.5,k2=0.5,k3=0.5,A=1)

# fitting
fitval=nls.lm(par=myparms,fn=ssq)

一旦我运行它，就会出现这样的错误

Error in chol.default(object$hessian) : 
  the leading minor of order 1 is not positive definite

Answer 1

你的代码的问题如下：

在代码行 cinit=c(A=myparms[4],B=0,C=0) 中，您为 A 提供了 myparms[4] 的值和 myparms[4] 的名称。让我们看看：

myparms=c(k1=0.5,k2=0.5,k3=0.5,A=1)
cinit=c(A=myparms[4],B=0,C=0)
print(cinit)
A.A   B   C 
1   0   0 

要解决这个问题，你可以这样做：

myparms=c(k1=0.5,k2=0.5,k3=0.5,A=1)
cinit=c(A=unname(myparms[4]),B=0,C=0)
print(cinit)
A   B   C 
1   0   0 

或者这个：

myparms=c(k1=0.5,k2=0.5,k3=0.5,1)
cinit=c(A=unname(myparms[4]),B=0,C=0)
print(cinit)
A   B   C 
1   0   0 

那么你的代码就可以工作了！

最好的问候， J_F