<!DOCTYPE HTML PUBLIC "-//W3C//DTD HTML 4.01 Transitional//EN">
<html>
<head>
  <meta content="text/html; charset=ISO-8859-1"
 http-equiv="Content-Type">
  <title></title>
</head>
<body bgcolor="#ffffff" text="#000000">
On 06/28/2010 07:36 AM, Ralph Kube wrote:
<blockquote cite="mid:4C289756.4070406@googlemail.com" type="cite">
  <pre wrap="">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 
<a class="moz-txt-link-freetext" href="http://www.tau.ac.il/~kineret/amit/scipy_tutorial/">http://www.tau.ac.il/~kineret/amit/scipy_tutorial/</a> in my code and it 
runs successfully.

I would be glad for any suggestion.


Cheers, Ralph
_______________________________________________
SciPy-User mailing list
<a class="moz-txt-link-abbreviated" href="mailto:SciPy-User@scipy.org">SciPy-User@scipy.org</a>
<a class="moz-txt-link-freetext" href="http://mail.scipy.org/mailman/listinfo/scipy-user">http://mail.scipy.org/mailman/listinfo/scipy-user</a>
  </pre>
</blockquote>
There are other optimization functions available such as those in the
'optimize' subpackage (
<meta http-equiv="content-type" content="text/html; charset=ISO-8859-1">
<a href="http://docs.scipy.org/scipy/docs/scipy.optimize/">http://docs.scipy.org/scipy/docs/scipy.optimize/)</a>
that may be more suited to this problem.<br>
<br>
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.<br>
<br>
Bruce<br>
</body>
</html>