[SciPy-user] use of Jacobian function in leastsq

[JJ]

Hello:
This is a question based on previous discussions by
Brendan Simons and Travis (RE: optimize.leastsq ->
fjac and the covariance redux).  I am wondering why
use of the Jacobian function in the leastsq function
does not seem to work.  There is probably a very
simple answer.  I am trying to use a Jacobian function
for a different problem and am using this as a test. 
The code is as follows, and you can see that C is not
equal to C2:

from scipy import *
from scipy.optimize import leastsq

x = array([0., 1., 2., 3., 4., 5.]).astype('d')
coeffs =  [4., 3., 5., 2.]
yErrs = array([0.1, 0.2, -0.1, 0.05, 0,
-.02]).astype('d')
m = len(x)
n = len(coeffs)

def evalY(p, x):
    """a simple cubic"""
    return x**3 * p[3] + x**2 * p[2] + x * p[1] + p[0]

def J(p, y, x):
    """the jacobian of a cubic (not actually
    a function of y, p since resid is linear in p)
    """
    print '\n in J'
    result = zeros([m,n], 'd')
    result[:, 0] = ones(m, 'd')
    result[:, 1] = x
    result[:, 2] = x**2
    result[:, 3] = x**3
    print result
    return result

def resid(p, y, x):
    return y - evalY(p, x)

y_true = evalY(coeffs, x)
y_meas = y_true + yErrs

p0 = [3., 2., 4., 1.]

pMin, C, infodict, ier, mesg = leastsq(resid, p0,
        args=(y_meas, x),  full_output=True)

pMin2, C2, infodict2, ier2, mesg2 = leastsq(resid, p0,
Dfun=J,
        args=(y_meas, x),  full_output=True)

#check against know algorithms for computing C
JAtPMin = J(pMin, None, x)
Q_true, R_true = linalg.qr(JAtPMin)
C_true = linalg.inv(dot(transpose(JAtPMin), JAtPMin))
C_true2 = linalg.inv(dot(transpose(R_true), R_true))

print 'pMin'
print pMin
print 'pMin2'
print pMin2
print '\nC'
print C
print '\nC2'
print C2
print '\nC_true'
print C_true
print '\nC_true2'
print C_true2

Thanks, JJ


 
____________________________________________________________________________________
Looking for earth-friendly autos? 
Browse Top Cars by "Green Rating" at Yahoo! Autos' Green Center.
http://autos.yahoo.com/green_center/