通过Scipy&Numpy使用Python将数据拟合到ODE系统

我通过Scipy& amp将我的MATLAB代码翻译成Python时遇到了一些麻烦. NumPy的.我坚持如何找到我的ODE系统的最佳参数值(k0和k1),以适应我的十个观察数据点.我目前对k0和k1有一个初步猜测.在MATLAB中,我可以使用一种叫做“fminsearch”的东西,它是一个接受ODE系统,观察数据点和ODE系统初始值的函数.然后,它将计算一对新的参数k0和k1,它们将适合观察到的数据.我已经包含了我的代码,看看你是否可以帮助我实现某种“fminsearch”来找到适合我数据的最佳参数值k0和k1.我想将任何代码添加到我的lsqtest.py文件中.

我有三个.py文件 – ode.py,lsq.py和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,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] #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,data)
    print resid

最佳答案

以下对我有用：

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,10)
y_data = np.array([0.000,0.309])


def f(y,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[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,teta)
    # calculate ode solution,retuen values for each entry of "x"
    r = integrate.odeint(f2,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] #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),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()