Sunday, July 7, 2013

Estimating ODE's parameters

In a previous post I used R to solve Ordinary Differential Equations (ODE). Since I'm still more familiar with Matlab to work with ODE, I decided to go one step further and learn how to use R to estimate parameters in ODE.

In this short tutorial, I will use the ODE presented here which quantifies the salt concentration in a well-stirred tank:

\[ \begin{align} &\frac{dx}{dt}=a\frac{150-x(t)}{200}\\[10pt] \end{align} \]

With \[ x(0) = 20 \]

The analytic solution is:

\[ \begin{align} &x(t)=150-130e^{-at/200} \end{align} \]

Step 1

Generate data using the analytic solution and add some random noise. Note that the ODE uses one parameter that has been fixed at \( a = 0.75 \).

t = seq(1, 800, by = 10)
a = 0.75  ## Fixed parameter used to simulate data.

## Simulate data and add noise.
conc.observed = 150 - 130 * exp(-a * t/200) + rnorm(length(t), sd = 5)

## Plot
plot(t, conc.observed, pch = 21, bg = "gray", ylim = c(20, 180), xlab = "Time", 
    ylab = "Salt (kg)", axes = F)
box(bty = "l")
axis(1, tck = -0.02)
axis(2, tck = -0.02, las = 2)

legend("bottomright", legend = c("Observed data"), pch = 21, pt.bg = "gray", 
    bty = "n")

plot of chunk chunk_8a

Step 2

Write a function that will be used to solve the ODE. While we there, solve it with fixed parameter \( a = 0.75 \).

library(deSolve)

salttank = function(t, state, parameters) {
    with(as.list(c(state, parameters)), {



        # rate of change
        dX = a * ((150 - X)/200)


        # return the rate of change
        list(c(dX))
    })
}


## Solve the ODE. Not necessary at this point.

## Define the initial value for the state variable
state = c(X = 20)

## Time range.
times = seq(1, 800, by = 10)

## Parameters
parameters = c(a = 0.75)

conc.modeled = ode(y = state, times = times, func = salttank, parms = parameters, 
    method = "rk4")

Plot the results.

plot(t, conc.observed, pch = 21, bg = "gray", ylim = c(20, 180), xlab = "Time", 
    ylab = "Salt (kg)", axes = F)
box(bty = "l")
axis(1, tck = -0.02)
axis(2, tck = -0.02, las = 2)
lines(conc.modeled[, "time"], conc.modeled[, "X"], col = "red")

legend("bottomright", legend = c("Observed data", "True solution (a = 0.75)"), 
    lty = c(NA, 1), pch = c(21, NA), col = c("black", "red"), pt.bg = "gray", 
    bty = "n")

plot of chunk chunk_8c

Step 3

This is where the magic happens. To fit parameters, I will use nls.lm() from the minpack.lm package. In its simplest form, the function uses par which is a list of guessed parameters and fn a function used to minimize the sum square of the residuals using the Levenberg-Marquardt algorithm.

Now we have to write a function (here I named it ssq) that will use ODE parameters as input (in this case only \( a \)) and produces a residuals vector as output.

ssq = function(params) {

    ## Range and initial condition as before.
    times = seq(1, 800, by = 10)
    state = c(X = 20)

    ## Resolve the ODE.
    out = ode(y = state, times = times, func = salttank, parms = params, method = "rk4")

    ## modeled - observed
    ssq = out[, "X"] - conc.observed
}

Step 4

Finally, we can estimate the parameter \( a \).


library(minpack.lm)

## Start with a guess a = 1
params.guessed = c(a = 1)
params.fitted = nls.lm(par = params.guessed, fn = ssq)

summary(params.fitted)
## 
## Parameters:
##   Estimate Std. Error t value Pr(>|t|)    
## a   0.7438     0.0113      66   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 4.93 on 79 degrees of freedom
## Number of iterations to termination: 5 
## Reason for termination: Relative error in the sum of squares is at most `ftol'.

We see from the summary of params.fitted that the estimated parameter is \( a = 0.7438 \) which is obviously very close to \( a=0.75 \) used to produce the observed data.

Step 5

Finally, simulate the data using the estimated parameter and plot the results.

## Simulate using fitted parameters
params = coef(params.fitted)

conc.modeled3 = ode(y = state, times = times, func = salttank, parms = params, 
    method = "rk4")
plot(t, conc.observed, pch = 21, bg = "gray", ylim = c(20, 180), xlab = "Time", 
    ylab = "Salt (kg)", axes = F)
box(bty = "l")
axis(1, tck = -0.02)
axis(2, tck = -0.02, las = 2)
lines(conc.modeled[, "time"], conc.modeled[, "X"], col = "red")
lines(conc.modeled3[, "time"], conc.modeled3[, "X"], col = "blue")

legend("bottomright", legend = c("Observed data", "True solution (a = 0.75)", 
    paste("Simulated with fitted parameters (a = ", sprintf("%2.2f", params), 
        ")", sep = "")), lty = c(NA, 1, 1), pch = c(21, NA, NA), col = c("black", 
    "red", "blue"), pt.bg = "gray", bty = "n")

plot of chunk chunk_8g

4 comments:

  1. Super,
    Je pensais que tu aurais fait cette modélisation pour la croissance d'une population de levure de bière! Peut être pour une prochaine.

    Sinon quand est-ce que tu t'attaque à PARAFAC?

    Bravo pour le travail de ton blog!
    Thomas

    ReplyDelete
  2. Merci Thomas!

    L'idée de la bière serait très intéressante. D'ailleurs, cela ferait un excellent sujet pour un cours à la maîtrise avec les étudiants. Suivre la cinétique des levures dans une expérience de fermentation.

    Bon courage,
    Phil

    ReplyDelete
  3. Many thanks for making the truthful effort to explain this. I feel very strong about it and would like to read more. If you
    can, as you find out more in depth knowledge, would you mind posting more posts similar to this one with more information.
    Qassim University
    .,

    ReplyDelete
    Replies
    1. Thank you for your comment. I have some posts coming soon with more detailed examples.

      Delete