fitting an ODE with python leastsq gives a cast error when initial conditions is passed as parameter

I have a set of data that I am trying to fit to an ODE model using scipy's leastsq function. My ODE has parameters beta and gamma, so that it looks for example like this:

# dS/dt = -betaSI
# dI/dt = betaSI - gammaI
# dR/dt = gammaI
# with y0 = y(t=0) = (S(0),I(0),R(0))

The idea is to find beta and gamma so that the numerical integration of my system of ODE's best approximates the data. I am able to do this just fine using leastsq if I know all the points in my initial condition y0.

Now, I am trying to do the same thing but to pass now one of the entries of y0 as an extra parameter. Here is where the Python and me stop communicating... I did a function so that now the first entry of the parameters that I pass to leastsq is the initial condition of my variable R. I get the following message:

*Traceback (most recent call last):
  File "/Users/Laura/Dropbox/SHIV/shivmodels/test.py", line 73, in <module>
    p1,success = optimize.leastsq(errfunc, initguess, args=(simpleSIR,[y0[0]],[Tx],[mydata]))
  File "/Library/Frameworks/Python.framework/Versions/7.2/lib/python2.7/site-packages/scipy/optimize/minpack.py", line 283, in leastsq
    gtol, maxfev, epsfcn, factor, diag)
TypeError: array cannot be safely cast to required type*

Here is my code. It is a little more involved that what it needs to be for this example because in reality I want to fit another ode with 7 parameters and want to fit to several data sets at once. But I wanted to post here something simpler... Any help will be very very much appreciated! Thank you very much!

import numpy as np
from matplotlib import pyplot as plt
from scipy import optimize
from scipy.integrate import odeint

#define the time span for the ODE integration:
Tx = np.arange(0,50,1)
num_points = len(Tx)

#define a simple ODE to fit:
def simpleSIR(y,t,params):
    dydt0 = -params[0]*y[0]*y[1]
    dydt1 = params[0]*y[0]*y[1] - params[1]*y[1]
    dydt2 = params[1]*y[1]
    dydt = [dydt0,dydt1,dydt2]
    return dydt


#generate noisy data:
y0 = [1000.,1.,0.]
beta = 12*0.06/1000.0
gamma = 0.25
myparam = [beta,gamma]
sir = odeint(simpleSIR, y0, Tx, (myparam,))

mydata0 = sir[:,0] + 0.05*(-1)**(np.random.randint(num_points,size=num_points))*sir[:,0]
mydata1 = sir[:,1] + 0.05*(-1)**(np.random.randint(num_points,size=num_points))*sir[:,1]
mydata2 = sir[:,2] + 0.05*(-1)**(np.random.randint(num_points,size=num_points))*sir[:,2]
mydata = np.array([mydata0,mydata1,mydata2]).transpose()


#define a function that will run the ode and fit it, the reason I am doing this
#is because I will use several ODE's to see which one fits the data the best.
def fitfunc(myfun,y0,Tx,params):
    myfit = odeint(myfun, y0, Tx, args=(params,))
    return myfit

#define a function that will measure the error between the fit and the real data:
def errfunc(params,myfun,y0,Tx,y):
    """
    INPUTS:
    params are the parameters for the ODE
    myfun is the function to be integrated by odeint
    y0 vector of initial conditions, so that y(t0) = y0
    Tx is the vector over which integration occurs, since I have several data sets and each 
    one has its own vector of time points, Tx is a list of arrays.
    y is the data, it is a list of arrays since I want to fit to multiple data sets at once
    """
    res = []
    for i in range(len(y)):
        V0 = params[0][i]
        myparams = params[1:]
        initCond = np.zeros([3,])
        initCond[:2] = y0[i]
        initCond[2] = V0
        myfit = fitfunc(myfun,initCond,Tx[i],myparams)
        res.append(myfit[:,0] - y[i][:,0])
        res.append(myfit[:,1] - y[i][:,1])
        res.append(myfit[1:,2] - y[i][1:,2])
    #end for
    all_residuals = np.hstack(res).ravel()
    return all_residuals
#end errfunc


#example of the problem:

V0 = [0]
params = [V0,beta,gamma]  
y0 = [1000,1]

#this is just to test that my errfunc does work well.
errfunc(params,simpleSIR,[y0],[Tx],[mydata])
initguess = [V0,0.5,0.5]

p1,success = optimize.leastsq(errfunc, initguess, args=(simpleSIR,[y0[0]],[Tx],[mydata]))

Solution

The problem is with the variable initguess. The function optimize.leastsq has the following call signature:

http://docs.scipy.org/doc/scipy/reference/generated/scipy.optimize.leastsq.html

It's second argument, x0, has to be an array. Your list

initguess = [v0,0.5,0.5]

won't be converted to an array because v0 is a list instead of an int or float. So you get an error when you try to convert initguess from a list to an array in the leastsq function.

I would adjust the variable params from

def errfunc(params,myfun,y0,Tx,y):

so that it is a 1-D array. Make the first few entries the values of v0 then append beta and gamma to that.