Fitting data to system of ODEs using Python via Scipy & Numpy

Question

I am having some trouble translating my MATLAB code into Python via Scipy & Numpy. I am stuck on how to find optimal parameter values (k0 and k1) for my system of ODEs to fit to my ten observed data points. I currently have an initial guess for k0 and k1. In MATLAB, I can using something called 'fminsearch' which is a function that takes the system of ODEs, the observed data points, and the initial values of the system of ODEs. It will then calculate a new pair of parameters k0 and k1 that will fit the observed data. I have included my code to see if you can help me implement some kind of 'fminsearch' to find the optimal parameter values k0 and k1 that will fit my data. I want to add whatever code to do this to my lsqtest.py file.

I have three .py files - ode.py, lsq.py, and lsqtest.py

ode.py:

def f(y, t, k): 
return (-k[0]*y[0],
      k[0]*y[0]-k[1]*y[1],
      k[1]*y[1])

lsq.py:

import pylab as py
import numpy as np
from scipy import integrate
from scipy import optimize
import ode
def lsq(teta,y0,data):
    #INPUT teta, the unknowns k0,k1
    # data, observed 
    # y0 initial values needed by the ODE
    #OUTPUT lsq value

    t = np.linspace(0,9,10)
    y_obs = data #data points
    k = [0,0]
    k[0] = teta[0]
    k[1] = teta[1]

    #call the ODE solver to get the states:
    r = integrate.odeint(ode.f,y0,t,args=(k,))


    #the ODE system in ode.py
    #at each row (time point), y_cal has
    #the values of the components [A,B,C]
    y_cal = r[:,1] #separate the measured B
    #compute the expression to be minimized:
    return sum((y_obs-y_cal)**2)

lsqtest.py:

import pylab as py
import numpy as np
from scipy import integrate
from scipy import optimize
import lsq


if __name__ == '__main__':

    teta = [0.2,0.3] #guess for parameter values k0 and k1
    y0 = [1,0,0] #initial conditions for system
    y = [0.000,0.416,0.489,0.595,0.506,0.493,0.458,0.394,0.335,0.309] #observed data points
    data = y
    resid = lsq.lsq(teta,y0,data)
    print resid

Solution 1

For these kind of fitting tasks you could use the package lmfit. The outcome of the fit would look like this; as you can see, the data are reproduced very well:

For now, I fixed the initial concentrations, you could also set them as variables if you like (just remove the vary=False in the code below). The parameters you obtain are:

[[Variables]]
    x10:   5 (fixed)
    x20:   0 (fixed)
    x30:   0 (fixed)
    k0:    0.12183301 +/- 0.005909 (4.85%) (init= 0.2)
    k1:    0.77583946 +/- 0.026639 (3.43%) (init= 0.3)
[[Correlations]] (unreported correlations are <  0.100)
    C(k0, k1)                    =  0.809 

The code that reproduces the plot looks like this (some explanation can be found in the inline comments):

import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import odeint
from lmfit import minimize, Parameters, Parameter, report_fit
from scipy.integrate import odeint


def f(y, t, paras):
    """
    Your system of differential equations
    """

    x1 = y[0]
    x2 = y[1]
    x3 = y[2]

    try:
        k0 = paras['k0'].value
        k1 = paras['k1'].value

    except KeyError:
        k0, k1 = paras
    # the model equations
    f0 = -k0 * x1
    f1 = k0 * x1 - k1 * x2
    f2 = k1 * x2
    return [f0, f1, f2]


def g(t, x0, paras):
    """
    Solution to the ODE x'(t) = f(t,x,k) with initial condition x(0) = x0
    """
    x = odeint(f, x0, t, args=(paras,))
    return x


def residual(paras, t, data):

    """
    compute the residual between actual data and fitted data
    """

    x0 = paras['x10'].value, paras['x20'].value, paras['x30'].value
    model = g(t, x0, paras)

    # you only have data for one of your variables
    x2_model = model[:, 1]
    return (x2_model - data).ravel()


# initial conditions
x10 = 5.
x20 = 0
x30 = 0
y0 = [x10, x20, x30]

# measured data
t_measured = np.linspace(0, 9, 10)
x2_measured = np.array([0.000, 0.416, 0.489, 0.595, 0.506, 0.493, 0.458, 0.394, 0.335, 0.309])

plt.figure()
plt.scatter(t_measured, x2_measured, marker='o', color='b', label='measured data', s=75)

# set parameters including bounds; you can also fix parameters (use vary=False)
params = Parameters()
params.add('x10', value=x10, vary=False)
params.add('x20', value=x20, vary=False)
params.add('x30', value=x30, vary=False)
params.add('k0', value=0.2, min=0.0001, max=2.)
params.add('k1', value=0.3, min=0.0001, max=2.)

# fit model
result = minimize(residual, params, args=(t_measured, x2_measured), method='leastsq')  # leastsq nelder
# check results of the fit
data_fitted = g(np.linspace(0., 9., 100), y0, result.params)

# plot fitted data
plt.plot(np.linspace(0., 9., 100), data_fitted[:, 1], '-', linewidth=2, color='red', label='fitted data')
plt.legend()
plt.xlim([0, max(t_measured)])
plt.ylim([0, 1.1 * max(data_fitted[:, 1])])
# display fitted statistics
report_fit(result)

plt.show()

If you have data for additional variables, you can simply update the function residual.

Solution 2

The following worked for me:

import pylab as pp
import numpy as np
from scipy import integrate, interpolate
from scipy import optimize

##initialize the data
x_data = np.linspace(0,9,10)
y_data = np.array([0.000,0.416,0.489,0.595,0.506,0.493,0.458,0.394,0.335,0.309])


def f(y, t, k): 
    """define the ODE system in terms of 
        dependent variable y,
        independent variable t, and
        optinal parmaeters, in this case a single variable k """
    return (-k[0]*y[0],
          k[0]*y[0]-k[1]*y[1],
          k[1]*y[1])

def my_ls_func(x,teta):
    """definition of function for LS fit
        x gives evaluation points,
        teta is an array of parameters to be varied for fit"""
    # create an alias to f which passes the optional params    
    f2 = lambda y,t: f(y, t, teta)
    # calculate ode solution, retuen values for each entry of "x"
    r = integrate.odeint(f2,y0,x)
    #in this case, we only need one of the dependent variable values
    return r[:,1]

def f_resid(p):
    """ function to pass to optimize.leastsq
        The routine will square and sum the values returned by 
        this function""" 
    return y_data-my_ls_func(x_data,p)
#solve the system - the solution is in variable c
guess = [0.2,0.3] #initial guess for params
y0 = [1,0,0] #inital conditions for ODEs
(c,kvg) = optimize.leastsq(f_resid, guess) #get params

print "parameter values are ",c

# fit ODE results to interpolating spline just for fun
xeval=np.linspace(min(x_data), max(x_data),30) 
gls = interpolate.UnivariateSpline(xeval, my_ls_func(xeval,c), k=3, s=0)

#pick a few more points for a very smooth curve, then plot 
#   data and curve fit
xeval=np.linspace(min(x_data), max(x_data),200)
#Plot of the data as red dots and fit as blue line
pp.plot(x_data, y_data,'.r',xeval,gls(xeval),'-b')
pp.xlabel('xlabel',{"fontsize":16})
pp.ylabel("ylabel",{"fontsize":16})
pp.legend(('data','fit'),loc=0)
pp.show()

Solution 3

    # cleaned up a bit to get my head around it - thanks for sharing 
    import pylab as pp
    import numpy as np
    from scipy import integrate, optimize

    class Parameterize_ODE():
        def __init__(self):
            self.X = np.linspace(0,9,10)
            self.y = np.array([0.000,0.416,0.489,0.595,0.506,0.493,0.458,0.394,0.335,0.309])
            self.y0 = [1,0,0] # inital conditions ODEs
        def ode(self, y, X, p):
            return (-p[0]*y[0],
                     p[0]*y[0]-p[1]*y[1],
                               p[1]*y[1])
        def model(self, X, p):
            return integrate.odeint(self.ode, self.y0, X, args=(p,))
        def f_resid(self, p):
            return self.y - self.model(self.X, p)[:,1]
        def optim(self, p_quess):
            return optimize.leastsq(self.f_resid, p_guess) # fit params

    po = Parameterize_ODE(); p_guess = [0.2, 0.3] 
    c, kvg = po.optim(p_guess)

    # --- show ---
    print "parameter values are ", c, kvg
    x = np.linspace(min(po.X), max(po.X), 2000)
    pp.plot(po.X, po.y,'.r',x, po.model(x, c)[:,1],'-b')
    pp.xlabel('X',{"fontsize":16}); pp.ylabel("y",{"fontsize":16}); pp.legend(('data','fit'),loc=0); pp.show()

Comments

[Zack]

Not sure if minimize would work for my case. Would you care to provide some code to show how? I figured something like the lsq function would be more for what I want to do, which is find optimal parameter values to match my data.

[tacaswell]

Welcome to SO. Your answer would be improved if you added a bit more commentary.

[Blackbam]

Please explain your answers.

[Andres Valdez]

Very interesting answer. This answer integrates the ode. And compares the outputs... could it be also used for MCMC (Markov Chain Monte Carlo)???

[Cleb]

@AndresValdez: Thank you! Quite an old answer so I should check whether that is still the best way of doing it in lmfit... Regarding your question, I do not know, probably best to check the documentation.

[Cleb]

@MNewville: Can you help Andres?

[Andres Valdez]

in a MCMC scheme I would need to solve the same ODE many times. "testing all the possible combinations of parameters"..... How efficient is this integration method?....

[Andres Valdez]

I mean the residual is defined in terms of the ode solution because the data is the measured solution at some time steps...... Mcmc will do this many times..... I wonder how could this work if having nonlinear odes

[Andres Valdez]

Apparently works very good. Found this post people.duke.edu/~ccc14/sta-663/CalibratingODEs.html this lmfit package looks very promising