使用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