현재이 튜토리얼 (http://www.r-bloggers.com/learning-r-parameter-fitting-for-models-involving-differential-equations/)에 따라 Levenberg-Marquardt 루틴 (nls.lm)을 pkg-minpack.lm
으로 사용하여 ODE 기능 응답을 맞추려고합니다. 예에서 nls.lm이있는 ODE 모델에 맞는 매개 변수 및 초기 조건
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
SSQ = 함수 잘못 입력 (myparms) cinit의 =의 앞에 # C (A = myparms [4], B = 0, C = 0)이어야 삭제됨 – Hong