Ajuste de datos al sistema de ODE utilizando Python a través de Scipy & Numpy

Question

Tengo problemas para traducir mi código MATLAB a Python a través de Scipy & Numpy. Estoy atascado en cómo encontrar valores de parámetros óptimos (k0 y k1) para que mi sistema de ODE se ajuste a mis diez puntos de datos observados. Actualmente tengo una conjetura inicial para k0 y k1. En MATLAB, puedo usar algo llamado 'fminsearch' que es una función que toma el sistema de ODEs, los puntos de datos observados y los valores iniciales del sistema de ODEs. Luego calculará un nuevo par de parámetros k0 y k1 que se ajustarán a los datos observados. He incluido mi código para ver si puede ayudarme a implementar algún tipo de 'fminsearch' para encontrar los valores de parámetros óptimos k0 y k1 que se ajusten a mis datos. Quiero agregar cualquier código para hacer esto a mi archivo lsqtest.py.

Tengo tres archivos .py - ode.py, lsq.py y 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 

Answer 1

Mire el scipy.optimizemodule. La función minimize es bastante similar a fminsearch, y creo que ambas usan básicamente un algoritmo simplex para la optimización.

Answer 2

A continuación trabajó para mí:

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() 

Answer 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() 

Answer 4

Para este tipo de tareas de ajuste que podría utilizar el paquete lmfit. El resultado del ajuste se vería así; como se puede ver, los datos se reproducen muy bien:

Por ahora, fija las concentraciones iniciales, también se les podría fijar como variables si se desea (acaba de quitar la vary=False en el código de abajo) . Los parámetros que obtenga son:

[[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 

El código que reproduce la trama se parece a esto (alguna explicación puede encontrarse en los comentarios en línea):

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() 

Si tiene datos para las variables adicionales, simplemente puede actualizar la función residual.