Temat: optymalizacja ukladu równan rózniczkowych

Dzien dobry,

Próbuje optymalizowac uklad równan rózniczkowych ze znalezieniem wartosci poczatkowych i funkcjami zmiennych (cE i cP) zapisanych w data.frame. Próbowalam wykorzystac ten kod na innym systemie równan i kod dzialal, chociaz generowal warning message :

In r[2] = (1/k1) * 2.42 * ((1 + 0.26 * cP)/(1 + 0.004 * ... :
number of items to replace is not a multiple of replacement length.

Mam niekompletne dane (równania sa dla pL i L), które oznaczaja koncentracje substancji, funkcje zewnetrzne oznaczaja rózny poziom kazdego dnia (cE i cP dane, z których tylko dane dla stezenia cP byly oznaczane kazdego dnia). Chce oznaczyc wartosci najbardziej wrazliwych parametrów, pozostale wartosci parametrów pochodza z literatury (mój próbny kod "nie dzialal", gdy próbowalam oznaczyc wieksza liczbe parametrów, niz liczba równan, wylaczajac oznaczane wartosci poczatkowe).

Czy moglabym prosic o wskazówki, gdzie moge znalezc przyklady, wyjasnienia jak "naprawic" mój kod, przyklady jak optymalizowac podobne uklady r-n rózniczkowych?

Nie moge znalezc pomocnego wyjasnienia w internecie, probowalam szukac wskazowek, jak usunac error messages.

Dziekuje z góry za jakakolwiek podpowiedz.

[code]tutaj jest kod[/ #set working directory
setwd("~/R/wkspace")
#load libraries
library(ggplot2)
library(reshape2)
library(deSolve)
library(minpack.lm)
time=c(22,23,24,46,47,48)
cE=c(15.92,24.01,25.29,15.92,24.01,25.29)
cP=c(0.3,0.14,0.29,0.3,0.14,0.29)
cL=c(6.13,3.91,38.4,6.13,3.91,38.4)
#CE and CP values for all cycle
#CE=(3.55,2.8,6.65,9.95,9.34,12.27 ...missing data ...12.43,15.92,24.01,25.29)
#at time (0-5 and 20-23) &
#cP=(0.18,0.28,2.01,8.88,9.34,20.44,37.76,30.72,30.53,37.88,51.66,50.02,43.73,
#48.7,64.78,82.93,51.24,44.7,59.94,27.09,0.59,0.3,0.14,0.29) at time (0-23)
df<-data.frame(time,cE,cP,cL)
df
names(df)=c("time","cE","cP","cL")

#rate function
rxnrate=function(t,c,parms){

#rate constant passed through a list called
#k1=parms$k1
#k2=parms$k2
#k3=parms$k3
#k4=parms$k4
#k5=parms$k5
#k6=parms$k6
#k7=parms$k7
k1=parms$k1
k2=parms$k2

#k10=parms$k10
#c is the concentration of species
#derivatives dc/dt are computed below
r=rep(0,length(c))
r[1]=(500+(4500*cE^8)/(200^8+cE^8))/(1+cP/12.2)-2.42*((1+0.26*cP)/(1+0.004*cE))*c["pLH"]; #dRP_LH/dt
r[2]=(1/k1)*2.42*((1+0.26*cP)/(1+0.004*cE))*c["pLH"]-k2*c["L"] #dL/dt
return(list(r))
}


ssq=function(myparms){
#initial concentration
cinit=c(pLH=myparms[3],L=myparms[4])
cinit=c(pLH=unname(myparms[3]),L=unname(myparms[4]))
print(cinit)
#time points for which conc is reported
#include the points where data is available
t=c(seq(0,46,2),df$time)
t=sort(unique(t))
#parameters from the parameters estimation
#k1=500
#k2=4500 #myparms[2]
#k3=200 #myparms[3]
#k4=2.42 #myparms[4]
#k5=0.26 #myparms[5]
#k6=12.2 #myparms[6]
#k7=0.004 #myparms[7]
k8=#myparms[8]
k9= #myparms[9]
#k10=myparms[10]
#k11=myparms[11]
#k12=myparms[12]
#solve ODE for a given set of parameters
out=ode(y=cinit,times=t,func=rxnrate,parms=list(k1=k1,k2=k2))

#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 parameters
myparms=c(k1=55,k2=24,pLH=19,LH=3.22)


#fitting
fitval=nls.lm(par=myparms,fn=ssq)
#summary of fit
summary(fitval)
#estimated parameter
parest=as.list(coef(fitval))

dof=3*nrow(df)-2
dof
#mean error
ms=sqrt(deviance(fitval)/dof)
#variance Covariance Matrix
S=vcov(fitval)
S
# plot of predicted vs experimental data
#simulated predicted profile at estimated params
cinit=c(pLH=19,L=3.22)
t=seq(0,5,0.2)
parms=as.list(parest)
out=ode(y=cinit,times=t,func=rxnrate,parms=parms)
outdf=data.frame(out)
names(outdf)=c("time","ca_pred","cb_pred","cc_pred")

# Overlay predicted profile with experimental data
tmppred=melt(outdf,id.var=c("time"),variable.name="species",value.name="conc")
tmpexp=melt(df,id.var=c("time"),variable.name="species",value.name="conc")
p=ggplot(data=tmppred,aes(x=time,y=conc,color=species,linetype=species))+geom_line()
p=p+geom_line(data=tmpexp,aes(x=time,y=conc,color=species,linetype=species))
p=p+geom_point(data=tmpexp,aes(x=time,y=conc,color=species))
p=p+scale_linetype_manual(values=c(0,1,0,1,0,1))
p=p+scale_color_manual(values=rep(c("red","blue","green"),each=2))+theme_bw()
print(p)


]

Pozdrawiam,
Malgorzata