Hi Doreen, If you send me your code which you tried using PyDSTool, I will happily advise you on how to fix it. I'm sure it will be possible for me to fix it! -Rob On Mon, Mar 31, 2008 at 4:09 PM, Doreen Mbabazi <doreen@aims.ac.za> wrote:

Hi,

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,))

Doreen

Gabriel Gellner

...

...

I have as my fitting function a system of differential equations(3) and I am supposed to get the optimal parameters to use by fitting data to this system. I tried to use scipy.optimize.leastsq. The difficulty is with the function that calculates the difference between the data and the values from the fitting function. The result of the fitting function is a list with three values and yet I only need one of the values to be subtracted from the data value. Could you not just subtract the index of the returned value? Say it is the first value you need then just use value[0].

If not can you give a small example, it will make it easier to help.

Gabriel _______________________________________________ SciPy-user mailing list SciPy-user@scipy.org http://projects.scipy.org/mailman/listinfo/scipy-user

_______________________________________________ SciPy-user mailing list SciPy-user@scipy.org http://projects.scipy.org/mailman/listinfo/scipy-user

-- Robert H. Clewley, Ph. D. Assistant Professor Department of Mathematics and Statistics Georgia State University 720 COE, 30 Pryor St Atlanta, GA 30303, USA tel: 404-413-6420 fax: 404-651-2246 http://www.mathstat.gsu.edu/~matrhc http://brainsbehavior.gsu.edu/

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

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

Hi, I am still working at my code and I have not yet been able to get optimal parameters. The problem being that actually my fitting function which is a system of odes is evaluated at only one set of parameters(the initial) thus the final results I get is the initial set of parameters. I need to find a way to be able to find integrate the system of odes at each iteration as the set of parameters change. In other words, Is there a way to get the parameters at each iteration while using the leastsq function? If that is possible how can I integrate my system of odes using these sets of parameters at each iteration? Some help please. Code is below. from scipy import * from scipy.integrate import odeint from scipy.optimize import leastsq import scipy.io.array_import import Gnuplot,Gnuplot.funcutils gp=Gnuplot.Gnuplot(persist=1) # , debug=1 def residuals(p,V,t): err = V - S(t,p) return err #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 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] pname = (['lamda','d','k','delta','pi','c']) lamda_0 = 0.1; d_0 = 0.01; k_0 = 0.60e-3; delta_0 = 0.35; pi_0 = 800; c_0 = 3.0 p = array([lamda_0 , d_0, k_0, delta_0,pi_0,c_0]) 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,)) def S(t,p): v = y[:,2] return v pbest = leastsq(residuals, p, args=(V,t),maxfev=2000) print pbest[0] Doreen. Rob Clewley

Hi Doreen,

...

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

Minpack expects an array type to be returned. You're returning a python list. So try this:

def residuals(p,y,t): return array([V - S(t_val,p) for t_val in t])

Here I've also removed the unnecessary i index variable -- Python lets you iterate directly over the contents of the list t.

-Rob

-- Robert H. Clewley, Ph. D. Assistant Professor Department of Mathematics and Statistics Georgia State University 720 COE, 30 Pryor St Atlanta, GA 30303, USA

tel: 404-413-6420 fax: 404-651-2246 http://www.mathstat.gsu.edu/~matrhc http://brainsbehavior.gsu.edu/ _______________________________________________ SciPy-user mailing list SciPy-user@scipy.org http://projects.scipy.org/mailman/listinfo/scipy-user

Wed, 02 Apr 2008 11:22:34 +0200, Doreen Mbabazi wrote:

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.

I think you should try to construct your program piece-by-piece, testing each part one-by-one (eg. printing/displaying some intermediate results to see whether they are correct), ie. 1. Does the function f(y,t,p) function as expected? 2. Does the function S(t, p) function as expected? 3. Does the function residuals(p, y, t) function as expected? 4. Does the whole shebang work? I think in this case you would have noticed problems in phases 2 and 3 even before 4. [clip]

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

I'm not sure whether the indentation was messed up in the mail, but if the code reads like this, it probably won't work (y_list will always contain only a single value). Anyway, you'll likely be better off using indexing instead of constructing a list: def S(t, p): ys = odeint(f, initial_y, t, args =(p,)) return ys[:,2]

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))]

Also this looks a bit funny: you create a list of len(t) elements that are all the same. Also, each element is an 1-d array containing the values V[:] - S(t[0],p). Anyway, you probably meant (no need to use the loop, subtraction is defined also for arrays): def residuals(p, y, t): return V - S(t,p) In general, if you find yourself writing xrange(len(z)) in numpy code, this often indicates that you should think a second time what you are trying to do: there may be an easier and more efficient way.

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.

I think this says that the 'residuals' function didn't return an array (1- d?) of floats, which it indeed doesn't appear to do. -- Pauli Virtanen