[SciPy-User] scipy.optimize.leastsq question

Bruce Southey bsouthey@gmail....
Mon Jun 28 08:44:32 CDT 2010


On 06/28/2010 07:36 AM, Ralph Kube wrote:
> Hello people,
> I am having a problem using the leastsq routine. My goal is to
> determine three parameters r_i, r_s and ppw so that the residuals
> to a model function a(r_i, r_s, ppw) to a measurement are minimal.
> When I call the leastsq routine with a good guess of starting values, it
> iterates 6 times without changing the vales of the initial parameters
> and then exits without an error.
> The function a is very complicated and expensive to evaluate. Some
> evaluation is done by using the subprocess module of python. Can this
> pose a problem for the leastsq routine?
>
>
> This is in the main routine:
>
> import numpy as N
>
> for t_idx, t in enumerate(time_var):
>
> 	r_i = 300.
> 	r_s = 1.0
> 	ppw=1e-6
> 	sza = 70.
> 	wl = N.arange(300., 3001., 1.)		
>
> 	albedo_true = compute_albedo(r_i, r_s, ppw, sza, wl)
> 	# This emulates the measurement data
> 	albedo_meas = albedo_true + 0.01*N.random.randn(len(wl))
> 		
> 	print 'Optimizing albedo'
> 	p0 = [2.*r_i, 1.4*r_s, 4.*ppw]
> 	plsq2 = leastsq(albedo_residual, p0, args=(albedo_meas, sza,
> wl))
> 	print '... done: ', plsq2[0][0], plsq2[0][1], plsq2[0][2]
> 	albedo_model = compute_albedo(plsq2[0][0], plsq2[0][1], plsq2[0][2],
> sza, wl)
>
> The residual function:
> def albedo_residual(p, y, sza, wvl):
> 	r_i, r_s, ppw = p
> 	albedo = compute_albedo(r_i, r_s, ppw, sza, wvl)
> 	err = albedo - y
> 	print 'Albedo for  r_i = %4.0f, r_s = %4.2f, ppw = %3.2e \
> 		Residual squared: %5f' % (r_i, r_s, ppw, N.sum(err**2))
>
> 	return err
>
> The definition of the function a(r_i, r_s, ppw)
> def compute_albedo(radius_ice, radius_soot, ppw, sza, wvl):
>
> The output is:
> Optimizing albedo
> Albedo for r_i =  600, r_s = 1.40, ppw = 4.00e-06 		Residual squared:
> 0.973819
> Albedo for r_i =  600, r_s = 1.40, ppw = 4.00e-06 		Residual squared:
> 0.973819
> Albedo for r_i =  600, r_s = 1.40, ppw = 4.00e-06 		Residual squared:
> 0.973819
> Albedo for r_i =  600, r_s = 1.40, ppw = 4.00e-06 		Residual squared:
> 0.973819
> Albedo for r_i =  600, r_s = 1.40, ppw = 4.00e-06 		Residual squared:
> 0.973819
> Albedo for r_i =  600, r_s = 1.40, ppw = 4.00e-06 		Residual squared:
> 0.973819
> ... done:  600.0 1.4 4e-06
>
> To check for errors, I implemented the example code from
> http://www.tau.ac.il/~kineret/amit/scipy_tutorial/ in my code and it
> runs successfully.
>
> I would be glad for any suggestion.
>
>
> Cheers, Ralph
> _______________________________________________
> SciPy-User mailing list
> SciPy-User@scipy.org
> http://mail.scipy.org/mailman/listinfo/scipy-user
>    
There are other optimization functions available such as those in the 
'optimize' subpackage ( 
http://docs.scipy.org/scipy/docs/scipy.optimize/) 
<http://docs.scipy.org/scipy/docs/scipy.optimize/> that may be more 
suited to this problem.

You probably have a scaling issue because your 'r_i' parameter is huge 
compared to your 'ppw' parameter (300 vs 0.000001). This is really 
really important if you model is nonlinear. So please try to standardize 
your values so that the parameters have similar magnitude - even just 
division/multiplication by some power of 10 can make a huge difference. 
If these parameters are so different or you need 'leastsq' then you 
probably should try either grid searching or fixing one or two 
parameters at a time. This will at least give you an idea on the 
possible values.

Bruce
-------------- next part --------------
An HTML attachment was scrubbed...
URL: http://mail.scipy.org/pipermail/scipy-user/attachments/20100628/57e6377e/attachment-0001.html 


More information about the SciPy-User mailing list