Lin Shao wrote:
Try this:
import numpy as N import scipy.optimize as O
## Create my data to fit: x=N.arange(-10,10,dtype=N.float32) y=(x-1.234)**2+3.456 ## I want to fit y=p0*(x-p1)^2+p2
## Now I define my objective def obj_func(params, xx, yy, mode='col'): return params[0] * (xx-params[1])**2 + params[2] - yy ## 'mode' is dummy here, because I want to use it for my Jocobian
## Now define my Jacobian def Jacobian(params,xx,yy,mode='col'): J = N.empty((len(params),xx.size)) J[0] = (xx-params[1])**2 J[1] = -2*params[0]*(xx-params[1])
shouldn't that be -2*params[0]*(xx-params[1])*params[1] ?
J[2] = 1 if mode=='col': return J elif mode=='row: return J.transpose() else: raise ValueError, "Unkown mode %s in Jacobian()" % mode ## Hopefully I did my calculus correctly
## First, try without Jocobian given b=O.leastsq(obj_func, (1.,1.,2.), args=(x, y, 'col')) print b[0] ## the result is [ 1. 1. 2.] ## obviously a failure, but no warning is returned
Not really a failure - this is numerics. The step length for finding the numeric derivative is too small (epsfcn keyword arg). Try with epsfcn = 1e-10. Have look at the return message of leastsq, which is b[3] in case full_output=1.
## Second, try Jacobian with col_deriv=1 b=O.leastsq(obj_func, (1.,1.,2.), args=(x, y, 'col'), Dfun=Jacobian, full_output=1, col_deriv=1) ## full_output has to be 1 because of a seg fault bug I reported earlier
does not segfault here. unix, python2.5, scipy 0.5.3.dev3062, numpy 1.0.2.dev3484
print b[0] ## the result is [ 1.00000004 1.23400008 3.45599962] ## Good Job!
## Last, try Jacobian with col_deriv=0 b=O.leastsq(obj_func, (1.,1.,2.), args=(x, y, 'row'), Dfun=Jacobian, full_output=1, col_deriv=0) print b[0] ## the result is [ 1.02507777 1.222815 2.2651552 ] ## completely different result and wrong
I've no idea why this fails, though. Christian