[SciPy-user] leastsq not converging with tutorial example

[Kael Fischer]

Matt et al:

Thanks for your reply. I'm still in the leastsq positive control
stage, so for now lets just stick to stuff based on the tutorial and
another example (http://linuxgazette.net/115/andreasen.html). 
Incidentally, my own application is simple 3-parameter first-order
decay (i.e. A = (A0-B)*exp(-k/t) + B), and I'm sure if I can get the
examples working, leastsq will be adequate.  Although I am not sure
how to report the quality of the fit, or how to estimate the
uncertainty of the fit parameters.

Here 2 scripts that I hope the mailer doesn't mangle.  Using either
(one is a sine wave; the other is a double exponential), I get back
the initial parameters and the "The cosine of the angle between
func(x) and any column of the\n  Jacobian is at most 0.000000 in
absolute value" status message.  Anybody know what is going on.

# Based on Scipy leastsq section of the tutorial
from scipy import *
from scipy.optimize import leastsq
x = arange(0,6e-2,6e-2/30)
A,k,theta = 10, 1.0/3e-2, pi/6
y_true = A*sin(2*pi*k*x+theta)
y_meas = y_true + 2*randn(len(x))

def residuals(p, y, x):
    A,k,theta = p
    err = y-A*sin(2*pi*k*x+theta)
    return err

def peval(x, p):
    return p[0]*sin(2*pi*p[1]*x+p[2])

p0 = [8, 1/2.3e-2, pi/3]
print array(p0)

plsq = leastsq(residuals, p0, args=(y_meas, x),full_output=True )
print plsq

==== or another examle ====

# From kinfit.py see link above
# fits to a double exponential
#
# data file avalilable at http://linuxgazette.net/115/misc/andreasen/tgdata.dat
#
from scipy import *
from scipy.optimize import leastsq
import scipy.io.array_import
from scipy import gplt


def residuals(p, y, x):
        err = y-peval(x,p)
        return err

def peval(x, p):
        return p[0]*(1-exp(-(p[2]*x)**p[4])) + p[1]*(1-exp(-(p[3]*(x))**p[5] ))

filename='tgdata.dat'
data = scipy.io.array_import.read_array(filename)

y = data[:,1]
x = data[:,0]

A1_0=5
A2_0=3
k1_0=0.5
k2_0=0.04
n1_0=2
n2_0=1
pname = (['A1','A2','k1','k2','n1','n2'])
p0 = array([A1_0 , A2_0, k1_0, k2_0,n1_0,n2_0])
plsq = leastsq(residuals, p0, args=(y, x), maxfev=2000,full_output=True )
print plsq

=== End Python Code ===

I just rebuilt Scipy by hand and made sure all the correct libraries
are there...  It has the same problem (initially I was using the
FreeBSD port to build locally).

That said - there were problems with cephs/protos.h. The following
complex functions were definitions were defined elsewhere (?) and I
comment them out:

//extern void cexp ( cmplx *z, cmplx *w );
//extern void csin ( cmplx *z, cmplx *w );
//extern void ccos ( cmplx *z, cmplx *w );
//extern void ctan ( cmplx *z, cmplx *w );
//extern void casin ( cmplx *z, cmplx *w );
//extern void cacos ( cmplx *z, cmplx *w );
//extern void catan ( cmplx *z, cmplx *w );
//extern void csqrt ( cmplx *z, cmplx *w );

Rgds,
Kael


On 10/25/05, Matt Fishburn <fishburn at mit.edu> wrote:
> Kael Fischer wrote:
> > Hi:
> >
> > When running the leastsq example form the tutorial, the data are not
> > fit and the starting parameters are returned.  Looking a little
> > deeper....
>
> I've used the leastsq algorithm to fit to general degree polynomials and
> an erf function before - see code at:
>
> http://web.mit.edu/fishburn/www/021/funcFit.py
>
> I haven't tried a cosine fit, though - do you have more in-depth code,
> or did you just grab the code from the tutorial?  I've had pretty good
> luck with leastsq.
>
> -Matt Fishburn
>


--
Kael Fischer, Ph.D
DeRisi Lab - Univ. Of California San Francisco
415-514-4320