我通过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()