[R] Optimize parameters of ODE Problem which is solved numeric

[Ravi Varadhan]

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:

> PID
      sato      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.