Michael S.
2010-Oct-24 00:05 UTC
[R] Optimize parameters of ODE Problem which is solved numeric
Hi, I have a data-matrix:> PIDsato hrs fim health 214 3 4.376430 6.582958 5 193 6 4.361825 3.138525 6 8441 6 4.205771 3.835886 7 7525 6 4.284489 3.245139 6 6806 7 4.168926 2.821833 7 5682 7 1.788707 1.212653 7 5225 6 1.651463 1.436980 7 4845 6 1.692710 1.267359 4 4552 5 1.686448 1.220539 6 4282 6 1.579868 1.086296 6 75441 6 2.978580 1.338100 7 I want to solve the following system of ode (ord. differential equations) numerically (f.e. with euler) dm1/dt <- a*m1+b*m2+d*m3+e*m4 dm2/dt <- a*m1+b*m2+d*m3+e*m4 dm3/dt <- a*m1+b*m2+d*m3+e*m4 dm4/dt <- a*m1+b*m2+d*m3+e*m4 with following initial values: m1<- PID$sato[1] m2<- PID$hrs[1] m3<- PID$fim[1] m4<- PID$health[1] a,b,d,e are free coeffient. The parameters a,b,d,e are not fix. My goal is to find the optimal parameters dependent on the error of the numerical solution compared to the registered values of PID. I would like to use instead of fix parameters, intervals. And then, through a minimizing function (minimal error) evaluate the optimal parameters. In order to have a first try, I picked fix parameters. a=0.1 b=0.22 d= -0.3 e= -0.4 #parameter vector parameter <- c(a,b,d,e) ### Initial Values of the ODE AWPs <- function(PID){ m1<- PID$sato[1] m2<- PID$hrs[1] m3<- PID$fim[1] m4<- PID$health[1] return(c(m1,m2,m3,m4)) } AWP <- AWPs(PID) ## before using ode from Package deSolve SWB <- function (t, AWP, parameter) { ### it didn't save the names of the variables, thats why this strange ### lines (but inimportant) m1<-AWP[1] m2<-AWP[2] m3<-AWP[3] m4<-AWP[4] ### with(as.list(c(AWP,parameter)),{ # Rate of change dm1 <- a*m1+b*m2+d*m3+e*m4 dm2 <- a*m1+b*m2+d*m3+e*m4 dm3 <- a*m1+b*m2+d*m3+e*m4 dm4 <- a*m1+b*m2+d*m3+e*m4 return(list(c(dm1, dm2, dm3, dm4))) }) } ### Time of data times <- seq(1996,2006, by=1) require(deSolve) ## instead of ode use also euler out <- ode(y = AWP, times = times, func = SWB, parms = parameter) ### Calculating the errors of the solution into a matrix. Taking sum is ### easy. R1 <- matrix(0,11,4) for (i in (1:11)){ for(j in (1:4)){ R1[i,j]<-out[i,j+1]-PID[i,j] } } sum(R1) What I tried next was, with the package sets, to use the class of interval. require(sets)> a1<-interval(-2,2)## Making 4-dimensional cuboit (carthesian product)> a4<-a^4## and trying: out <- ode(y = AWP, times = times, func = SWB, parms = a4) But ode can not process, as it looks like, an interval object, since it is a list and not numeric. But I couldn't manage to make a4 or a1 numeric with as.numeric(a4). The results of "out" don't make much sense. MY QUESTIONS: 1. Is it possible to use an interval instead of fix variables for ode or euler? 2. Reading all of it, is there maybe a easier way to find the optimal parameters for a,b,d,e 3. Does somebody have an idea which function would minimize something like: min{g(x)-f(x)|x in [-2,2]} as you see g(x) is a bit complicated. Thanks for your help in advance. Michael -- GMX DSL Doppel-Flat ab 19,99 €/mtl.! Jetzt auch mit gratis Notebook-Flat! http://portal.gmx.net/de/go/dsl
Ravi Varadhan
2010-Oct-25 14:00 UTC
[R] Optimize parameters of ODE Problem which is solved numeric
Hi Michael, You do not need a numerical solver for this. This is a linear system of ODEs and it admits closed form solutions. The solution is given as: Y(t) = c_1 * v_1 * exp(k_1 * t) + ... + c_4 * v_4 * exp(k_4 * t) where k_1, ..., k_4 are the eigenvalues (can be real or complex) and v_1, ..., v_4 are the eigenvectors (each is a vector of length 4) of the 4x4 coefficient matrix . You can get the constants c_1, ..., c_4 by solving for the initial condition. However, your specification of the coefficient matrix is incorrect. You can cannot have the same coefficients in all the equations. Ravi. -----Original Message----- From: r-help-bounces at r-project.org [mailto:r-help-bounces at r-project.org] On Behalf Of Michael S. Sent: Saturday, October 23, 2010 8:06 PM To: r-help at r-project.org Subject: [R] Optimize parameters of ODE Problem which is solved numeric Hi, I have a data-matrix:> PIDsato hrs fim health 214 3 4.376430 6.582958 5 193 6 4.361825 3.138525 6 8441 6 4.205771 3.835886 7 7525 6 4.284489 3.245139 6 6806 7 4.168926 2.821833 7 5682 7 1.788707 1.212653 7 5225 6 1.651463 1.436980 7 4845 6 1.692710 1.267359 4 4552 5 1.686448 1.220539 6 4282 6 1.579868 1.086296 6 75441 6 2.978580 1.338100 7 I want to solve the following system of ode (ord. differential equations) numerically (f.e. with euler) dm1/dt <- a*m1+b*m2+d*m3+e*m4 dm2/dt <- a*m1+b*m2+d*m3+e*m4 dm3/dt <- a*m1+b*m2+d*m3+e*m4 dm4/dt <- a*m1+b*m2+d*m3+e*m4 with following initial values: m1<- PID$sato[1] m2<- PID$hrs[1] m3<- PID$fim[1] m4<- PID$health[1] a,b,d,e are free coeffient. The parameters a,b,d,e are not fix. My goal is to find the optimal parameters dependent on the error of the numerical solution compared to the registered values of PID. I would like to use instead of fix parameters, intervals. And then, through a minimizing function (minimal error) evaluate the optimal parameters. In order to have a first try, I picked fix parameters. a=0.1 b=0.22 d= -0.3 e= -0.4 #parameter vector parameter <- c(a,b,d,e) ### Initial Values of the ODE AWPs <- function(PID){ m1<- PID$sato[1] m2<- PID$hrs[1] m3<- PID$fim[1] m4<- PID$health[1] return(c(m1,m2,m3,m4)) } AWP <- AWPs(PID) ## before using ode from Package deSolve SWB <- function (t, AWP, parameter) { ### it didn't save the names of the variables, thats why this strange ### lines (but inimportant) m1<-AWP[1] m2<-AWP[2] m3<-AWP[3] m4<-AWP[4] ### with(as.list(c(AWP,parameter)),{ # Rate of change dm1 <- a*m1+b*m2+d*m3+e*m4 dm2 <- a*m1+b*m2+d*m3+e*m4 dm3 <- a*m1+b*m2+d*m3+e*m4 dm4 <- a*m1+b*m2+d*m3+e*m4 return(list(c(dm1, dm2, dm3, dm4))) }) } ### Time of data times <- seq(1996,2006, by=1) require(deSolve) ## instead of ode use also euler out <- ode(y = AWP, times = times, func = SWB, parms = parameter) ### Calculating the errors of the solution into a matrix. Taking sum is ### easy. R1 <- matrix(0,11,4) for (i in (1:11)){ for(j in (1:4)){ R1[i,j]<-out[i,j+1]-PID[i,j] } } sum(R1) What I tried next was, with the package sets, to use the class of interval. require(sets)> a1<-interval(-2,2)## Making 4-dimensional cuboit (carthesian product)> a4<-a^4## and trying: out <- ode(y = AWP, times = times, func = SWB, parms = a4) But ode can not process, as it looks like, an interval object, since it is a list and not numeric. But I couldn't manage to make a4 or a1 numeric with as.numeric(a4). The results of "out" don't make much sense. MY QUESTIONS: 1. Is it possible to use an interval instead of fix variables for ode or euler? 2. Reading all of it, is there maybe a easier way to find the optimal parameters for a,b,d,e 3. Does somebody have an idea which function would minimize something like: min{g(x)-f(x)|x in [-2,2]} as you see g(x) is a bit complicated. Thanks for your help in advance. Michael -- GMX DSL Doppel-Flat ab 19,99 €/mtl.! Jetzt auch mit gratis Notebook-Flat! http://portal.gmx.net/de/go/dsl ______________________________________________ R-help at r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.