[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