[R-sig-dyn-mod] Error in r[1] <- : replacement has length zero and Error in chol.default(object$hessian) : the leading minor of order 1 is not positive definite

[Malgorzata Wieteska]

Hello,
I would like to ask for a help. I'm quite new to R, I try to find initial conditions while optimising system of differential equations. Unfortunately I get an error messages: Error in r[1] <- :    replacement has length zero and Error in chol.default(object$hessian) :    the leading minor of order 1 is not positive definite

Thanks in advance for any help. 
My code:
library(ggplot2)library(reshape2)library(deSolve)library(minpack.lm)time<-seq(0, 5, by =1)F=c(0.408,0.226,0.126,0.224,0.123,0.162)signal<-as.data.frame(list(time=time,FSH))input<-approxfun(signal,rule=2)cB=c(0,9491.0179,6779.799,1186.7998,0,0)
df<-data.frame(time,F,cSeF)dfnames(df)=c("time","F","cB")#plot datatmp=melt(df,id.vars=c("time"), variable.name="species",value.name="conc")ggplot(data=tmp,aes(x=time,y=conc,color=species))+geom_point(size=3)#rate functionrxnrate=function(t,c,parms){  #rate constant passed through a list called  a=parms$a  c2=parms$c2  k1=parms$k1  c3=parms$c3  p1=parms$p1  #k4=parms$k4  #c is the concentration of species  #derivatives dc/dt are computed below  FSH<-input(t)  r[1]<-a*F+k1*F*c["A"]- c2*c["A"]#dA/dt  r[2]<-c2*c["A"]- c3*c["B"]+p1*F*c["B"]^2 #dA/dt    res<-list(r)  list(res,signal=FSH)    return(list(r)) }# predicted concentration for a given parametercinit=c(A=452,B=0)t=df$timeparms=list(a=0.0107,k1=0.321,c2=0.749,c3=0.749,p1=0.214)out=ode(y=cinit,times=t,func=rxnrate,parms=list(a=a,k1=k1,c2=c2,c3=c3,p1=p1))head(out)
ssq=function(myparms){   #initial concentration  cinit=c(A=myparms[6],B=myparms[7])  cinit=c(A=unname(myparms[6]),B=unname(myparms[7]))  print(cinit)  #time points for which conc is reported  #include the points where data is available  t=c(seq(0,5,1),df$time)  t=sort(unique(t))  #parameters from the parameters estimation  a=myparms[1]  c1=myparms[2]  c2=myparms[3]  c3=myparms[4]  p1=myparms[5]  #k5=myparms[5]  #solve ODE for a given set of parameters  out=ode(y=cinit,time=t,func=rxnrate,parms=list(a=a,c1=c1,c2=c2,c3=c3,p1=p1))  #Filter data that contains time points  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(ssqres)}# parameter fitting using levenberg marquart#initial guess for parametersmyparms=c(a=0.0107,k1=0.321,c2=0.749,c3=0.749,p1=0.214,A=452,B=0)

#fittingfastpath=FALSEfitval=nls.lm(par=myparms,fn=ssq)#summary of fitsummary(fitval)#estimated parameterparest=as.list(coef(fitval))#degrees of freedom: # data points=#paramdof=2*nrow(df)-2dof#mean errorms=sqrt(deviance(fitval)/dof)#variance Covariance MatrixS=vcov(fitval)S# plot of predicted vs experimental data#simulated predicted profile at estimated paramscinit=c(myparms[6],myparms[7])t=seq(0,5,1)parms=as.list(parest)out=ode(y=cinit,times=t,func=rxnrate,parms=parms)outdf=data.frame(out)names(outdf)=c("time","A_pred","B_pred")
	[[alternative HTML version deleted]]