[SciPy-user] Estimation of parameters while fitting data

Anne Archibald peridot.faceted@gmail....
Mon Mar 31 16:51:42 CDT 2008

On 31/03/2008, Doreen Mbabazi <doreen@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,

More information about the SciPy-user mailing list