[SciPy-User] scipy.optimize.leastsq question
Ralph Kube
ralphkube@googlemail....
Mon Jun 28 07:36:38 CDT 2010
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
More information about the SciPy-User
mailing list