[SciPy-user] Estimation of parameters while fitting data
Doreen Mbabazi
doreen@aims.ac...
Wed Apr 2 04:22:34 CDT 2008
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@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@scipy.org
> http://projects.scipy.org/mailman/listinfo/scipy-user
>
More information about the SciPy-user
mailing list