title: Fitting data to system of ODEs using Python via Scipy & Numpy
body: <p>I am having some trouble translating my MATLAB code into Python via Scipy &amp; Numpy. I am stuck on how to find optimal parameter values (k0 and k1) for my system of ODEs to fit to my ten observed data points. I currently have an initial guess for k0 and k1. In MATLAB, I can using something called 'fminsearch' which is a function that takes the system of ODEs, the observed data points, and the initial values of the system of ODEs. It will then calculate a new pair of parameters k0 and k1 that will fit the observed data. I have included my code to see if you can help me implement some kind of 'fminsearch' to find the optimal parameter values k0 and k1 that will fit my data. I want to add whatever code to do this to my lsqtest.py file. </p>
  
  <p>I have three .py files - ode.py, lsq.py, and lsqtest.py </p>
  
  <p>ode.py: </p>
  
  <pre><code>def f(y, t, k): 
  return (-k[0]*y[0],
        k[0]*y[0]-k[1]*y[1],
        k[1]*y[1])
  </code></pre>
  
  <p>lsq.py: </p>
  
  <pre><code>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)
  </code></pre>
  
  <p>lsqtest.py: </p>
  
  <pre><code>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
  </code></pre>