[SciPy-user] Fitting with global parameters optimize.leastsq

[Steve Schmerler]

Iain Day wrote:
>>> I'm messing about with scipy.reshape at the moment and sort of getting 
>>> somewhere. The residuals is now a 1D array, as is t and y. However, when 
>>> I use leastsq I get the following error:
>>>
>>>   File "/sw/lib/python2.4/site-packages/scipy/optimize/minpack.py", line 
>>> 229, 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.
>> Could you post an example code producing this error?
>>
> 
> I hope the line wrapping doesn't hamper things too much. The model 
> function is working fine, and the simple fit of the model to one trace 
> works fine, its the extension to simutaneous fitting of several traces 
> thats got me stumped.
> 
> Thanks for your help,
> 
> 
> 
> import scipy, scipy.optimize, scipy.io
> 
> rawdata = scipy.io.read_array(open(datadir + filename + ".csv", "rb"), 
> separator=',')
> 
> data_points = rawdata[:, 0]
> raw_signals = rawdata[:, 1:]
> 
> # These are the guesses. I want to use the same fitted value to
> # B across all traces, it should be a common parameter.
> A = [40.0, 40.0, 40.0, 40.0, 40.0, 40.0]
> B = 0.3
> C = [45.0, 45.0, 45.0, 3.0, 3.0, 15.0]
> 
> p0 = scipy.r_[A, B, C]
> fix0 = [x, y, z]
> 
> expdata = scipy.reshape(raw_signals, (1, len(A)*len(data_points)))
> 
> p1, success = scipy.optimize.leastsq(residuals, p0, args=(fix0, 
> data_points, expdata))
> 
> 
> def residuals(varpars, fixpars, points, expdata):
>      ntraces = (len(varpars) - 1) / 2
>      A = varpars[:ntraces]
>      B = varpars[ntraces]
>      C = varpars[ntraces+1:]
> 
>      err = scipy.zeros((len(points), ntraces), dtype=float)
>      raw_data = scipy.zeros((len(points), ntraces), dtype=float)
>      raw_data = scipy.reshape(expdata, (len(points), ntraces))
> 
>      for i in range(ntraces):
>          err[:, i] = raw_data[:, i] - model([A[i], B, C[i]], fixpars, 
> points)
> 
>      err = scipy.reshape(err, (1, ntraces*len(points)))
> 
>      return err
> 

What's x, y, z and def model(...) is missing :). I thought more of a 
copy-paste-and-run example (e.g. with rawdata = numpy.random.randn(100, 
7) where 100 == len(data_points)) that produces the error :-)

Anyway:

1) The error comes from the shape of err in your residuals(). leastsq 
expects an array of shape (<number>,), not (1, <number>) (i.e. 
len(err.shape) = 1, not 2.
You would have to do reshape(err, (ntraces*len(points),)).

I get this error when I replace

def residuals(x, t, y):
     return y - model(t, x)

by

def residuals(x, t, y):
     err =  y - model(t, x)
     err = scipy.reshape(err, (1, len(err)))
     print "### err.shape:", err.shape
     return err


2) You're trying to fit all data sets (t == data_points, y == 
raw_signals[:, i]) at once and I guess that won't work. You have to loop 
over A (or C or raw_signals[:, i]) building one data set and one initial 
guess (x0 = numpy.array([A[i], B, C[i]])) per loop.


-- 
cheers,
steve

Random number generation is the art of producing pure gibberish as 
quickly as possible.