nls.lmを使用したODEモデルのパラメータと初期条件

Question

ここでは、チュートリアル（http://www.r-bloggers.com/learning-r-parameter-fitting-for-models-involving-differential-equations/）に続いて、Levenberg-Marquardtルーチン（nls.lm）をpkg-minpack.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)あなたは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