[SciPy-user] scipy.optimize.leastsq col_deriv strangeness w/ data fit

[John Hunter]

I am using leastsq to do a best fit to a simple exponential function.
In my test script, I find that if I use col_deriv=0, I get a different
answer than I get if I do a col_deriv=1 with the err func returning
the transposed jacobian.  When I compare the true parameters with the
best fit parameters, the correct answer is with col_deriv=1 but not
with col_deriv=0.  Am I misusing this parameter or is something amiss?

In the data fitting problem, x is the independent variable with
parameters parsTrue.  I create the data values by assuming they are 

  data = func(pars) + noise

Here is the test script that works, to fit the data to 

where func(pars) is a*exp(alpha*x) + k


To get the best fit params, I minimize over pars

  data - func(pars)

------------------------------------------------------------
from scipy import exp, arange, zeros, Float, ones
from RandomArray import normal
from scipy.optimize import leastsq

parsTrue = 2.0, -.76, 0.1   # a, alpha, k
x = arange(0, 4, 0.001)

def func(pars):
    a, alpha, k = pars
    return a*exp(alpha*x) + k

def errfunc(pars):
    return data - func(pars)  #return the error


def deriv_errfunc(pars):
    'The Jacobian of the errfunc is -Jacobian of the func'

    a, alpha, k = pars
    J = zeros( (3,len(x)), Float)  # init the Jacobian only once
    #deriv with respect to a
    J[0,:] = -exp(alpha*x)

    #deriv with respect to alpha
    J[1,:] = -x*a*exp(alpha*x)

    #deriv with respect to k
    J[2,:] = -ones((len(x),), Float) # this is indep of k
    return J


# some pseudo data; add some noise
data = func(parsTrue) + normal(0.0, 0.1, x.shape)


guess = 1.0, -.4, 0.0   # the intial guess of the params

best, info, ier, mesg = leastsq(errfunc, guess,
                                full_output=1,
                                Dfun=deriv_errfunc,
                                col_deriv=1,
                                )

print 'true', parsTrue
print 'best', best

------------------------------------------------------------
This returns best fit params
true (2.0, -0.76000000000000001, 0.10000000000000001)
best [ 2.00487303 -0.74859602  0.09035577]



However, If I replace the deriv_errfunc with row derivs, manually
transposing the jacobian at creation time

def deriv_errfunc(pars):
    'The Jacobian of the errfunc is -Jacobian of the func'

    a, alpha, k = pars
    J = zeros( (len(x),3), Float)  # init the Jacobian only once

    # deriv with respect to a
    J[:,0] = -exp(alpha*x)

    # deriv with respect to alpha
    J[:,1] = -x*a*exp(alpha*x)

    # deriv with respect to k
    J[:,2] = -ones((len(x),), Float) 
    return J


and 

best, info, ier, mesg = leastsq(errfunc, guess,
                                full_output=1,
                                Dfun=deriv_errfunc,
                                col_deriv=0,
                                )
The fit is poor, eg, 

true (2.0, -0.76000000000000001, 0.10000000000000001)
best [ 1.13708025 -0.40303463  0.1669024 ]


With col_deriv=0 sometimes the solution is found and sometimes not
depending on the run, but it is consistently around the values
reported above (ie, wrong).  With col_deriv=1, though, the solution is
always found and is correct.


Suggestions?

John Hunter