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

Question

我在通过 Scipy 和 Numpy 将我的 MATLAB 代码转换为 Python 时遇到了一些问题。我坚持如何为我的 ODE 系统找到最佳参数值(k0 和 k1)以适合我的十个观察数据点。我目前对 k0 和 k1 有一个初步猜测。在 MATLAB 中，我可以使用称为“fminsearch”的函数，它是一个函数，它采用 ODE 系统、观察到的数据点和 ODE 系统的初始值。然后它将计算一对新的参数 k0 和 k1，以适应观察到的数据。我已经包含了我的代码，看看您是否可以帮助我实现某种“fminsearch”来找到适合我的数据的最佳参数值 k0 和 k1。我想将执行此操作的任何代码添加到我的 lsqtest.py 文件中。

我有三个 .py 文件 - ode.py、lsq.py 和 lsqtest.py

代码.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

对于这些类型的装配任务，您可以使用包 lmfit .适合的结果看起来像这样；如您所见，数据重现得很好:

目前，我固定了初始浓度，如果您愿意，您也可以将它们设置为变量(只需删除下面代码中的 vary=False)。你得到的参数是:

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

重现情节的代码如下所示(可以在内联注释中找到一些解释):

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

如果您有其他变量的数据，您可以简单地更新函数 residual。