[SciPy-User] Alternatives to scipy.optimize
Mon Feb 27 09:56:30 CST 2012
On Sun, Feb 26, 2012 at 10:14 AM, Matthew Newville
> Hi Erik, Josef,
> On Saturday, February 25, 2012 8:21:43 PM UTC-6, Erik Petigura wrote:
>> Thanks for getting back to me!
>> I'd like to minimize p1 and p2 together. Let me try to describe my
>> problem a little better:
>> I'm trying to fit an exoplanet transit light curve. My model is a box + a
>> polynomial trend.
>> The polynomial coefficients and the depth of the box are linear
>> parameters, so I want to
>> fit them using linear least squares. The center and width of the transit
>> are non-linear
>> so I need to fit them with an iterative approach like optimize.fmin.
>> Here's how I implemented it.
> I'm not sure I fully follow your model, but if I understand correctly,
> you're looking to find optimal parameters for something like
> model = linear_function(p1) + nonlinear_function(p2)
yes, I've read about this mostly in the context of semiparametric
versions when the nonlinear function does not have a parametric form.
> for sets of coefficients p1 and p2, each set having a few fitting variables,
> some of which may be related. Is there an instability that prevents you
> from just treating this as a single non-linear model?
I think p1 and p2 shouldn't have any cross restrictions or it will get
a bit more complicated.
The main reason for splitting it up is computational, I think. It's
quite common in econometrics to "concentrate out" parameters that have
an explicit solution, so we need the nonlinear optimization only for a
smaller parameter space.
> Another option might be to have the residual function for
> scipy.optimize.leastsq (or lmfit) call numpy.linalg.lstsq at each
> iteration. I would think that more fully explore the parameter space than
> first fitting nonlinear_function with scipy.optimize.fmin() then passing
> those best-fit parameters to numpy.linalg.lstsq(), but perhaps I'm not fully
> understanding the nature of the problem.
That's what both of us did, my version https://gist.github.com/1911544
>> There is a lot unpacking and repacking the parameter array as it gets
>> passed around
>> between functions. One option that might work would be to define
>> functions based on a
>> "parameter object". This parameter object could have attributes like
>> linear/non-linear. I found a more object oriented optimization module
>> However, it doesn't allow for linear fitting.
> Linear fitting could probably be added to lmfit, though I haven't looked
> into it. For this problem, I would pursue the idea of treating your
> fitting problem as a single model for non-linear least squares with
> optimize.leastsq or with lmfit. Perhaps I missing something about your
> model that makes this approach unusually challenging.
> Josef P wrote:
>> The easiest is to just write some helper functions to stack or unstack
>> the parameters, or set some to fixed. In statsmodels we use this in
>> some cases (as methods since our models are classes), also to
>> transform parameters.
>> Since often this affects groups of parameters, I don't know if the
>> lmfit approach would helps in this case.
> If many people who are writing their own model functions find themselves
> writing similar helper functions to stack and unstack parameters, "the
> easiest" here might not be "the best", and providing tools to do this
> stacking and unstacking might be worthwhile. Lmfit tries to do this.
>> (Personally, I like numpy arrays with masks or fancy indexing, which
>> is easy to understand. Ast manipulation scares me.)
> I don't understand how masks or fancy indexing would help here. How would
> that work?
---- <detour warning>
a few examples how we currently handle parameter restrictions in statsmodels
Skippers example: transform parameters in the loglikelihood to force
the parameter estimates of an ARMA process to produce a stationary
(stable) solution - transforms groups of parameters at once
not a clean example: fitting distributions with some frozen parameters
select parameters that are not frozen to use with fmin
x0 = np.array(x0)[np.isnan(frmask)]
expand the parameters again to include the frozen parameters inside
the loglikelihood function
theta = frmask.copy()
theta[np.isnan(frmask)] = thetash
It's not as user friendly as the version that got into
scipy.stats.distribution, (but it's developer friendly because I don't
have to stare at it for hours to spot a bug)
structural Vector Autoregression uses a similar mask pattern
In another sandbox model I build a nested dictionary to map the model
parameters to the reduced list that can be fed to
But we don't have user defined nonlinear constraints yet.
> FWIW, lmfit uses python's ast module only for algebraic constraints between
> parameters. That is,
> from lmfit import Parameter
> Parameter(name='a', value=10, vary=True)
> Parameter(name='b', expr='sqrt(a) + 1')
> will compile 'sqrt(a)+1' into its AST representation and evaluate that for
> the value of 'b' when needed. So lmfit doesn't so much manipulate the AST
> as interpret it. What is manipulated is the namespace, so that 'a' is
> interpreted as "look up the current value of Parameter 'a'" when the AST is
> evaluated. Again, this applies only for algebraic constraints on
It's a bit similar to a formula framework for specifying a statistical
model that is to be estimated (with lengthy discussion on the
I see the advantages but I haven't spent the weeks of time to figure
out what's behind the machinery that is required (especially given all
the other statistics and econometrics that is missing, and where I
only have to worry about how numpy, scipy and statsmodels behaves. And
I like Zen of Python #2)
> Having written fitting programs that support user-supplied algebraic
> constraints between parameters in Fortran77, I find interpreting python's
> AST to be remarkably simple and robust. I'm scared much more by statistical
> modeling of economic data ;)
different tastes and applications.
I'd rather think about the next 20 statistical tests I want to code,
than about the AST or how sympy translates into numpy.
> --Matt Newville
> SciPy-User mailing list
More information about the SciPy-User