2

我目前正在尝试在实验室实验中使用 minpack 中的 Levenberg-Marquardt 例程 (nls.lm) 来拟合功能响应，以解决消耗问题。作为一个例子，我一直在按照这里的教程（http://www.r-bloggers.com/learning-r-parameter-fitting-for-models-involving- ）在 minpack 中使用 levenberg-marquardt 例程（nls.lm）微分方程/）。

在示例中，他首先设置了一个我修改过的函数 rxnrate 来拟合数据，如下所示：

# rate function
rxnrate=function(t,c,parms){

# rate constant passed through a list called parms
a=parms$a
h=parms$h

# c is the concentration of species

# derivatives are computed below
r=rep(0,length(c))
r[1]=-c["B"]*a*c["A"]/(c["B"]+a*h*c["A"])#prey
r[2]=0#predator


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

}

我的问题是，我没有使用很长的时间序列，而是有许多具有多个起点的短时间序列（n=6）。用 nls.lm 函数单独拟合这些将导致相当无用的估计。我的低技术解决方案与罗杰斯分析方法产生了可比的结果，它是将它们全部排列并同时拟合它们，如下例所示。

# rate function
rxnrate=function(t,c,parms){

 # rate constant passed through a list called parms
 a=parms$a
 h=parms$h

 # c is the concentration of species

 # derivatives are computed below
 r=rep(0,length(c))
 r[1]=-c["B"]*a*c["A"]/(c["B"]+a*h*c["A"])#prey
 r[2]=0#predator
 r[3]=-c["D"]*a*c["C"]/(c["D"]+a*h*c["C"])#prey2
 r[4]=0#predator2
 r[5]=-c["F"]*a*c["E"]/(c["F"]+a*h*c["E"])#prey3
 r[6]=0#predator3
# and so on

 return(list(r))

}

这样做的问题是我很快就用完了字母，而且硬编码所有这些时间序列（总共超过 100 个）效率极低。

我的问题是因为配对方程都是相同的，是否有一个解决方案可以让我编写一次并让函数将其应用于所有后续配对时间序列。我也不确定这个解决方案是否会导致参数估计的任何数学问题，即使它似乎给出了与其他方法相当的结果。

这是一个小的工作示例

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

    #load population data
rate= structure(list(time = c(0, 0.5, 1, 1.5, 2, 2.5), a = c(6L, 5L, 
      3L, 4L, 3L, 3L), b = c(1L, 1L, 1L, 1L, 1L, 1L), c = c(6L, 3L, 
      3L, 4L, 2L, 3L), d = c(3L, 3L, 3L, 3L, 3L, 3L), e = c(6L, 6L, 
      4L, 2L, 3L, 3L), f = c(6L, 6L, 6L, 6L, 6L, 6L), g = c(12L, 8L, 
      8L, 8L, 8L, 7L), h = c(1L, 1L, 1L, 1L, 1L, 1L), i = c(12L, 11L, 
      7L, 6L, 3L, 4L), j = c(3L, 3L, 3L, 3L, 3L, 3L), k = c(24L, 12L, 
      11L, 15L, 8L, 7L), l = c(1L, 1L, 1L, 1L, 1L, 1L)), .Names = c("time", 
      "a", "b", "c", "d", "e", "f", "g", "h", "i", "j", "k", "l"), row.names = c(NA, 
      6L), class = "data.frame")


rxnrate=function(t,c,parms){

  # rate constant passed through a list called parms
  a=parms$a
  h=parms$h  
  m=parms$m

  # derivatives dc/dt are computed below
  r=rep(0,length(c))
  holling<-c["B"]*a*c["A"]/(c["B"]+a*h*c["A"]) 
  r[1]=-c["B"]*a*c["A"]/(c["B"]^m+a*h*c["A"]) #dN1/dt
  r[2]=0
  r[3]=-c["D"]*a*c["C"]/(c["D"]^m+a*h*c["C"]) #dN1/dt
  r[4]=0
  r[5]=-c["F"]*a*c["E"]/(c["F"]^m+a*h*c["E"]) #dN1/dt
  r[6]=0
  r[7]=-c["H"]*a*c["G"]/(c["H"]^m+a*h*c["G"]) #dN1/dt
  r[8]=0
  r[9]=-c["J"]*a*c["I"]/(c["J"]^m+a*h*c["I"]) #dN1/dt
  r[10]=0
  r[11]=-c["L"]*a*c["K"]/(c["L"]^m+a*h*c["K"]) #dN1/dt
  r[12]=0
  return(list(r))

}

ssq=function(parms){

  # inital concentration
  cinit=cinit
  # time points for which conc is reported
  # include the points where data is available
  t=c(seq(0,2.5,0.5),rate$time)
  t=sort(unique(t))
  # parameters from the parameter estimation routin
  a=parms[1]
  h=parms[2]
  m=parms[3]
  # solve ODE for a given set of parameters
  out=ode(y=cinit,times=t,func=rxnrate,parms=parms)

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

  return(ssqres)

}

 # parameter fitting using levenberg marquart algorithm
# initial guess for parameters
control=nls.lm.control(maxiter = 1000,ptol=0.000000000000000000000001,ftol=0.0000000000000000000001)
cinit=c(A=6,B=1,C=6,D=3,E=6,F=6,G=12,H=1,I=12,J=3,K=24,L=1)
parms=list(a=1,h=0.1,m=1)
fit=nls.lm(par=parms,fn=ssq,lower=c(rep(0,3)),upper=c(2,0.5,2),control=control)
4

0 回答 0