[SciPy-user] Estimation of parameters while fitting data

[Doreen Mbabazi]

Thanks,

What you proposed has worked and I have been able to get the values that I
want to use to obtain the residuals however when it comes to using the
leastsq, I get an error which I have tried to google but its giving me no
useful results. Below is the code and the error.

initial_y = [10,0,10e-6]
def f(y,t,p):
    y_dot = [0.,0.,0.]
    y_dot[0] = p[0] - p[1]*y[0] - p[2]*y[0]*y[2]
    y_dot[1] = p[2]*y[0]*y[2] - p[3]*y[1]
    y_dot[2] = p[4]*y[1] - p[5]*y[2]
    return y_dot
filename=('essay1.dat')  # data file is stored
data = scipy.io.array_import.read_array(filename) # data file is read.

t = data[:,0]
V = data[:,1]


lamda_0 = 0.16; d_0 = 0.016;  k_0 = 0.40e-3;  delta_0 = 0.35;   pi_0 =
900;   c_0 = 3
#p0 = array([lamda_0 , d_0, k_0, delta_0,pi_0,c_0])
p = array([lamda_0 , d_0, k_0, delta_0,pi_0,c_0])

def S(t, p):
    y_list = []
    ys = odeint(f, initial_y, t, args =(p,))
   for i in range(len(t)):
	y_V = ys[i][2]
	y_list.append(y_V)
	return y_list

y = odeint(f,initial_y,t,args=(p,))
#print S(t,p)

def residuals(p,y,t):
    return [V - S(t,p) for i in xrange(len(t))]

pbest = leastsq(residuals, p, args=(V,t))

The error
ValueError: setting an array element with a sequence.
Traceback (most recent call last):
  File "essay6.py", line 44, in <module>
    pbest = leastsq(residuals, p, args=(V,t))
  File "/usr/lib/python2.5/site-packages/scipy/optimize/minpack.py", line
266, in leastsq
    retval =
_minpack._lmdif(func,x0,args,full_output,ftol,xtol,gtol,maxfev,epsfcn,factor,diag)
minpack.error: Result from function call is not a proper array of floats.

Doreen.

Anne Archibald
> On 31/03/2008, Doreen Mbabazi <doreen at aims.ac.za> wrote:
>
>>  Thanks, I tried to do that(by taking err = V-f(y,t,p)[2]) while
>> defining
>>  the function residuals but the trouble is that actually f(y,t,p)
>>  calculates value of y at t0 so it cannot help me. What I want are the
>>  third values from y(y[i][2]). Below I have tried to do that but that
>> gives
>>  particular values of y so my parameters are not optimized.
>>
>>  def residuals(p, V, t):
>>      """The function is used to calculate the residuals
>>      """
>>      for i in range(len(t)):
>>          err = V-y[i][2]
>>          return err
>>
>>  #Function defined with y[0]=T,y[1]=T*,y[2] = V,lamda = p[0],d = p[1],
>>  k=p[2],delta=p[3], pi = p[4], c = p[5]
>>  initial_y = [10,0,10e-6] # initial conditions T(0)= 10cells , T*(0)=0,
>>  V(0)=10e-6
>>
>>  p is the list of parameters that are being estimated
>> (lamda,d,k,delta,pi,c)
>>  def f(y,t,p):
>>     y_dot = [0,0,0]
>>     y_dot[0] = p[0] - p[1]*y[0] - p[2]*y[0]*y[2]
>>     y_dot[1] = p[2]*y[0]*y[2] - p[3]*y[1]
>>     y_dot[2] = p[4]*y[1] - p[5]*y[2]
>>     return y_dot
>>
>>  y = odeint(f,initial_y,t,args=(p,))
>
> First of all, I'm not totally sure I have correctly understood your
> problem. Let me state it as I understand it:
>
> You have a collection of points, (x[i], y[i]). You also have a model y
> = S(x, P[j]) giving y as a function of x and some parameters P. This
> function is not given explicitly, instead you know that it is the
> solution to a certain differential equation, with initial values and
> parameters given by P. You want to find the values of P that minimize
> sum((y[i] - S(x[i],P))**2).
>
> Is that about right? (The differential equation is actually expressed
> as three coupled ODEs, but that's not really a problem.)
>
> The easiest-to-understand way to solve the problem is probably to
> start by writing a python function S that behaves like S does. Of
> course, it has to be computed by solving the ODE, which means we're
> going to have to solve the ODE a zillion times, but that's okay,
> that's what computers are for.
>
> def S(x, P):
>      ys = odeint(f, initial_y, [0,x], P)
>      return ys[1,0]
>
> Now check that this function looks vaguely right (perhaps by plotting
> it, or checking that the values that come out are sensible).
>
> Now you can do quite ordinary least-squares fitting:
>
> def residuals(P,x,y):
>     return [y[i] - S(x[i],P) for i in xrange(len(x))]
>
> Pbest = scipy.optimize.leastsq(residuals, Pguess, args=(x,y))
>
> This should work, and be understandable. But it is not very efficient,
> since for every set of parameters, we solve the ODE len(x) times. We
> can improve things by using the fact that odeint can return the
> integrated function evaluated at a list of places. So we'd modify S to
> accept a list of xs, and return the S values at all those places. This
> would even simplify our residuals function:
>
> def S(xs, P):
>      ys = odeint(f, initial_y, numpy.concatenate(([0],xs), P)
>      return ys[1:,0]
>
> def residuals(P,xs,ys):
>      return ys - S(xs, P)
>
> Is this the problem you were trying to solve? I suggest first getting
> used to how odeint and leastsq work first, then combining them. Their
> arguments can be weird, in particular the way odeint treats the
> initial x like the xs you want your ode integrated to.
>
> Good luck,
> Anne
> _______________________________________________
> SciPy-user mailing list
> SciPy-user at scipy.org
> http://projects.scipy.org/mailman/listinfo/scipy-user
>