3

I'm trying to fit an ODE model to some data and solve for the values of the parameters in the model.

I know there is a package called FME in R which is designed to solve this kind of problem. However, when I tried to write the code like the manual of this package, the program could not run with the following traceback information:

Error in lsoda(y, times, func, parms, ...) : illegal input detected before taking any integration steps - see written message

The code is the following:

x <- c(0.1257,0.2586,0.5091,0.7826,1.311,1.8636,2.7898,3.8773)
y <- c(11.3573,13.0092,15.1907,17.6093,19.7197,22.4207,24.3998,26.2158)

time <- 0:7

# Initial Values of the Parameters
parms <- c(r = 4, b11 = 1, b12 = 0.2, a111 = 0.5, a112 = 0.1, a122 = 0.1)

# Definition of the Derivative Functions
# Parameters in pars; Initial Values in y
derivs <- function(time, y, pars){
     with(as.list(c(pars, y)),{
         dx <- r + b11*x + b12*y - a111*x^2 - a122*y^2 - a112*x*y
         dy <- r + b12*x + b11*y - a122*x^2 - a111*y^2 - a112*x*y
         list(c(dx,dy))
     })
}

initial <- c(x = x[1], y = y[1])

data <- data.frame(time = time, x = x, y = y)

 # Cost Computation, the Object Function to be Minimized
 model_cost <- function(pars){
     out <- ode(y = initial, time = time, func = derivs, parms = pars)
     cost <- modCost(model = out, obs = data, x = "time")
     return(cost)
 }

 # Model Fitting
 model_fit <- modFit(f = model_cost, p = parms, lower = c(-Inf,rep(0,5)))

Is there anyone that knows how to use the FME package and fix the problem?

Ken.w
  • 31
  • 3

1 Answers1

2

Your code-syntax is right and it works until the last line.

you can check your code with

model_cost(parms)

This works fine and you can see with

model_cost(parms)$model

that your "initial guess" is far away from the observed data (compare "obs" and "mod"). Perhaps here is a failure so that the fitting procedure will not reach the observed data.

So much for the while ... I also checked different methods with parameter "methods = ..." but still does not work.

Best wishes, Johannes

edit: If you use:

model_fit <- modFit(f = model_cost, p = parms)

without any lower bounds, then you will get an result (even there are warnings), but then a112 is negative what you wanted to omit.

J_F
  • 9,956
  • 2
  • 31
  • 55
  • Maybe I should try different initial values for the parameters, or just reconsider the range of the parameters. Thanks for your reply :) – Ken.w Apr 28 '16 at 05:34
  • Another comment: not all optimizers support constraints natively and especially do not allow a mix (that includes `-Inf` or `Inf`). A possible approach is shown in this [recent post](https://stackoverflow.com/a/75908176/3677576). – tpetzoldt Apr 20 '23 at 05:00