[SciPy-user] use of Jacobian function in leastsq

lechtlr lechtlr@yahoo....
Sun Mar 18 11:13:36 CDT 2007


JJ:

If you define the derivatives of the Jacobian function in columns instead of rows, it'll work better. Try the attached script.

-Lex

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([n,m], '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 = [1., 1., 1., 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, col_deriv=1, args=(y_meas, x),  full_output=True)

#check against known algorithms for computing C
JAtPMin = J(pMin, None, x)
Q_true, R_true = linalg.qr(JAtPMin)
C_true = linalg.inv(dot(JAtPMin, transpose(JAtPMin)))
C_true2 = linalg.inv(dot(R_true, transpose(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



 
---------------------------------
Don't be flakey. Get Yahoo! Mail for Mobile and 
always stay connected to friends.
-------------- next part --------------
An HTML attachment was scrubbed...
URL: http://projects.scipy.org/pipermail/scipy-user/attachments/20070318/82ee5416/attachment-0001.html 


More information about the SciPy-user mailing list