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

[Michael S.]

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