सिप्पी और न्यूम्पी

Question

के माध्यम से पाइथन का उपयोग करते हुए ओडीई की प्रणाली में फिटिंग डेटा को मुझे सिस्टी & नम्पी के माध्यम से पाइथन में अपने MATLAB कोड का अनुवाद करने में कुछ परेशानी हो रही है। मैं अपने दस मनाए गए डेटा पॉइंट्स में फिट होने के लिए ओडीई की प्रणाली के लिए इष्टतम पैरामीटर मान (k0 और k1) कैसे ढूंढ सकता हूं, इस पर अटक गया हूं। मैं वर्तमान में k0 और k1 के लिए प्रारंभिक अनुमान है। MATLAB में, मैं 'fminsearch' नामक कुछ का उपयोग कर सकता हूं जो एक ऐसा कार्य है जो ओडीई, मनाए गए डेटा बिंदुओं और ओडीई की प्रणाली के प्रारंभिक मानों की प्रणाली लेता है। इसके बाद पैरामीटर k0 और k1 की एक नई जोड़ी की गणना की जाएगी जो मनाए गए डेटा को फिट करेगी। मैंने यह देखने के लिए अपना कोड शामिल किया है कि क्या आप मेरे डेटा को फिट करने वाले इष्टतम पैरामीटर मान k0 और k1 को खोजने के लिए किसी प्रकार के 'fminsearch' को लागू करने में मेरी सहायता कर सकते हैं। मैं अपने 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,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

को देखेंmodule। minimize फ़ंक्शन fminsearch के समान दिखता है, और मेरा मानना ​​है कि दोनों मूल रूप से अनुकूलन के लिए एक सरल एल्गोरिदम का उपयोग करते हैं।

Answer 2

निम्नलिखित मेरे लिए काम किया:

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

आप पैकेज lmfit इस्तेमाल कर सकते हैं फिटिंग कार्य इस तरह के लिए। फिट का नतीजा इस तरह दिखेगा; जैसा कि आप देख सकते हैं, डेटा बहुत अच्छी तरह से reproduced रहे हैं:

अभी के लिए, मैं प्रारंभिक सांद्रता तय, आप भी चर के रूप में उन्हें सेट करता है, तो आप की तरह (बस 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 

कोड है कि reproduces साजिश (कुछ स्पष्टीकरण इनलाइन टिप्पणी में पाया जा सकता) इस तरह दिखता है:

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 फ़ंक्शन को अपडेट कर सकते हैं।