2016-09-06 2 views
2

현재이 튜토리얼 (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 모델에 맞는 매개 변수 및 초기 조건

그는 제 I는 다음과 같이 변형 된 함수 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 
+2

SSQ = 함수 잘못 입력 (myparms) cinit의 =의 앞에 # C (A = myparms [4], B = 0, C = 0)이어야 삭제됨 – Hong

답변

2

문제처럼 나오는 한 다음입니다 : 코드 라인에서

다음은 내 코드입니다 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

+0

안녕하세요, J_F 님, 답변 해 주셔서 감사합니다. 이제 제대로 작동합니다. – Hong

+0

제 답변을 투표하고 최선의 답변으로 표시 할 수 있습니다. 고맙습니다. –

+0

또한 이것을 사용할 수 있습니다 :'myparms = c (k1 = 0.5, k2 = 0.5, k3 = 0.5, 1)'. 이것이 기반이 된 인용 된 예제는 명명 된 벡터가 아닌 목록을 사용했고, R을 목록으로 제공하는 최적화 루틴의 매개 변수는 훨씬 더 일반적이라는 점에 유의해야합니다. –

관련 문제