Fitting generic models with least-squares minimization (pwkit.lsqmdl)

Model data with least-squares fitting

This module provides tools for fitting models to data using least-squares optimization.

There are four basic approaches all offering a common programming interface:

• Generic Nonlinear Modeling

• One-dimensional Polynomial Modeling

• Modeling of a Single Scale Factor

• Modeling With Pluggable Components

ModelBase(data[, invsigma])	An abstract base class holding data and a model for least-squares fitting.
Parameter(owner, index)	Information about a parameter in a least-squares model.

class pwkit.lsqmdl.ModelBase(data, invsigma=None)

An abstract base class holding data and a model for least-squares fitting.

The models implemented in this module all derive from this class and so inherit the attributes and methods described below.

A Parameter data structure may be obtained by indexing this object with either the parameter’s numerical index or its name. I.e.:

m = Model(...).solve(...)
p = m['slope']
print(p.name, p.value, p.uncert, p.uval)

chisq = None

After fitting, the χ² of the fit.

covar = None

After fitting, the variance-covariance matrix representing the parameter uncertainties.

data = None

The data to be modeled; an n-dimensional Numpy array.

invsigma = None

Data weights: 1/σ for each data point.

make_frozen_func(params)

Return a data-generating model function frozen at the specified parameters.

As with the mfunc attribute, the resulting function may or may not take arguments depending on the particular kind of model being evaluated.

mdata = None

After fitting, the modeled data at the best parameters.

mfunc = None

After fitting, a callable function evaluating the model fixed at best params.

The resulting function may or may not take arguments depending on the particular kind of model being evaluated.

params = None

After fitting, a Numpy ndarray of solved model parameters.

plot(modelx, dlines=False, xmin=None, xmax=None, ymin=None, ymax=None, **kwargs)

Plot the data and model (requires omega).

This assumes that data is 1D and that mfunc takes one argument that should be treated as the X variable.

pnames = None

A list of textual names for the parameters.

print_soln()

Print information about the model solution.

puncerts = None

After fitting, a Numpy ndarray of 1σ uncertainties on the model parameters.

rchisq = None

After fitting, the reduced χ² of the fit, or None if there are no degrees of freedom.

resids = None

After fitting, the residuals: resids = data - mdata.

set_data(data, invsigma=None)

Set the data to be modeled.

Returns self.

show_corr()

Show the parameter correlation matrix with pwkit.ndshow_gtk3.

show_cov()

Show the parameter covariance matrix with pwkit.ndshow_gtk3.

class pwkit.lsqmdl.Parameter(owner, index)

Information about a parameter in a least-squares model.

These data may only be obtained after solving least-squares problem. These objects reference information from their parent objects, so changing the parent will alter the apparent contents of these objects.

property index

The parameter’s index in the Model’s arrays.

property name

The parameter’s name.

property uncert

The uncertainty in value.

property uval

Accesses value and uncert as a pwkit.msmt.Uval.

property value

The parameter’s value.

Generic Nonlinear Modeling

Model(simple_func, data[, invsigma, args])	Models data with a generic nonlinear optimizer
Parameter(owner, index)	Information about a parameter in a least-squares model.

class pwkit.lsqmdl.Model(simple_func, data, invsigma=None, args=())

Models data with a generic nonlinear optimizer

Basic usage is:

def func(p1, p2, x):
    simulated_data = p1 * x + p2
    return simulated_data

x = [1, 2, 3]
data = [10, 14, 15.8]
mdl = Model(func, data, args=(x,)).solve(guess).print_soln()

The Model constructor can take an optional argument invsigma after data; it specifies inverse sigmas, not inverse variances (the usual statistical weights), for the data points. Since most applications deal in sigmas, take care to write:

m = Model(func, data, 1. / uncerts) # right!

not:

m = Model(func, data, uncerts) # WRONG

If you have zero uncertainty on a measurement, you must wind a way to express that constraint without including that measurement as part of the data vector.

lm_prob = None

A pwkit.lmmin.Problem instance describing the problem to be solved.

After setting up the data-generating function, you can access this item to tune the solver.

make_frozen_func(params)

Returns a model function frozen to the specified parameter values.

Any remaining arguments are left free and must be provided when the function is called.

For this model, the returned function is the application of functools.partial() to the func property of this object.

set_func(func, pnames, args=())

Set the model function to use an efficient but tedious calling convention.

The function should obey the following convention:

def func(param_vec, *args):
    modeled_data = { do something using param_vec }
    return modeled_data

This function creates the pwkit.lmmin.Problem so that the caller can futz with it before calling solve(), if so desired.

Returns self.

set_simple_func(func, args=())

Set the model function to use a simple but somewhat inefficient calling convention.

The function should obey the following convention:

def func(param0, param1, ..., paramN, *args):
    modeled_data = { do something using the parameters }
    return modeled_data

Returns self.

solve(guess)

Solve for the parameters, using an initial guess.

This uses the Levenberg-Marquardt optimizer described in pwkit.lmmin.

Returns self.

One-dimensional Polynomial Modeling

class pwkit.lsqmdl.PolynomialModel(maxexponent, x, data, invsigma=None)

Least-squares polynomial fit.

Because this is a very specialized kind of problem, we don’t need an initial guess to solve, and we can use fast built-in numerical routines.

The output parameters are named “a0”, “a1”, … and are stored in that order in PolynomialModel.params[]. We have y = sum(x**i * a[i]), so “a2” = “params[2]” is the quadratic term, etc.

This model does not give uncertainties on the derived coefficients. The as_nonlinear() method can be use to get a Model instance with uncertainties.

Methods:

as_nonlinear - Return a (lmmin-based) Model equivalent to self.

as_nonlinear(params=None)

Return a Model equivalent to this object. The nonlinear solver is less efficient, but lets you freeze parameters, compute uncertainties, etc.

If the params argument is provided, solve() will be called on the returned object with those parameters. If it is None and this object has parameters in self.params, those will be use. Otherwise, solve() will not be called on the returned object.

make_frozen_func(params)

Return a data-generating model function frozen at the specified parameters.

As with the mfunc attribute, the resulting function may or may not take arguments depending on the particular kind of model being evaluated.

Modeling of a Single Scale Factor

class pwkit.lsqmdl.ScaleModel(x, data, invsigma=None)

Solve data = m * x for m.

make_frozen_func(params)

Return a data-generating model function frozen at the specified parameters.

As with the mfunc attribute, the resulting function may or may not take arguments depending on the particular kind of model being evaluated.

Modeling With Pluggable Components

ComposedModel(component, data[, invsigma])
ModelComponent([name])
AddConstantComponent([name])
AddValuesComponent(nvals[, name])	XXX terminology between this and AddConstant is mushy.
AddPolynomialComponent(maxexponent, x[, name])
SeriesComponent([components, name])	Apply a set of subcomponents in series, isolating each from the other.
MatMultComponent(k[, name])	Given a component yielding k**2 data points and k additional components, each yielding n data points.
ScaleComponent([subcomp, name])

class pwkit.lsqmdl.ComposedModel(component, data, invsigma=None)

debug_derivative(guess)

returns (explicit, auto)

make_frozen_func()

Return a data-generating model function frozen at the specified parameters.

As with the mfunc attribute, the resulting function may or may not take arguments depending on the particular kind of model being evaluated.

mfunc(*args)

class pwkit.lsqmdl.ModelComponent(name=None)

deriv(pars, jac)

Compute the Jacobian. jac[i] is d`mdata`/d`pars[i]`.

extract(pars, perr, cov)

Extract fit results into the object for ease of inspection.

finalize_setup()

If the component has subcomponents, this should set their name, setguess, setvalue, and setlimit properties. It should also set npar (on self) to the final value.

model(pars, mdata)

Modify mdata based on pars.

prep_params()

This should make any necessary calls to setvalue or setlimit, though in straightforward cases it should just be up to the user to do this. If the component has subcomponents, their prep_params functions should be called.

class pwkit.lsqmdl.AddConstantComponent(name=None)

deriv(pars, jac)

Compute the Jacobian. jac[i] is d`mdata`/d`pars[i]`.

extract(pars, perr, cov)

Extract fit results into the object for ease of inspection.

model(pars, mdata)

Modify mdata based on pars.

class pwkit.lsqmdl.AddValuesComponent(nvals, name=None)

XXX terminology between this and AddConstant is mushy.

deriv(pars, jac)

Compute the Jacobian. jac[i] is d`mdata`/d`pars[i]`.

extract(pars, perr, cov)

Extract fit results into the object for ease of inspection.

model(pars, mdata)

Modify mdata based on pars.

class pwkit.lsqmdl.AddPolynomialComponent(maxexponent, x, name=None)

deriv(pars, jac)

Compute the Jacobian. jac[i] is d`mdata`/d`pars[i]`.

extract(pars, perr, cov)

Extract fit results into the object for ease of inspection.

model(pars, mdata)

Modify mdata based on pars.

class pwkit.lsqmdl.SeriesComponent(components=(), name=None)

Apply a set of subcomponents in series, isolating each from the other. This is only valid if every subcomponent except the first is additive – otherwise, the Jacobian won’t be right.

add(component)

This helps, but direct manipulation of self.components should be supported.

deriv(pars, jac)

Compute the Jacobian. jac[i] is d`mdata`/d`pars[i]`.

extract(pars, perr, cov)

Extract fit results into the object for ease of inspection.

finalize_setup()

If the component has subcomponents, this should set their name, setguess, setvalue, and setlimit properties. It should also set npar (on self) to the final value.

model(pars, mdata)

Modify mdata based on pars.

prep_params()

This should make any necessary calls to setvalue or setlimit, though in straightforward cases it should just be up to the user to do this. If the component has subcomponents, their prep_params functions should be called.

class pwkit.lsqmdl.MatMultComponent(k, name=None)

Given a component yielding k**2 data points and k additional components, each yielding n data points. The result is [A]×[B], where A is the square matrix formed from the first component’s output, and B is the (k, n) matrix of stacked output from the final k components.

Parameters are ordered in same way as the components named above.

deriv(pars, jac)

Compute the Jacobian. jac[i] is d`mdata`/d`pars[i]`.

extract(pars, perr, cov)

Extract fit results into the object for ease of inspection.

finalize_setup()

If the component has subcomponents, this should set their name, setguess, setvalue, and setlimit properties. It should also set npar (on self) to the final value.

model(pars, mdata)

Modify mdata based on pars.

prep_params()

This should make any necessary calls to setvalue or setlimit, though in straightforward cases it should just be up to the user to do this. If the component has subcomponents, their prep_params functions should be called.

class pwkit.lsqmdl.ScaleComponent(subcomp=None, name=None)

deriv(pars, jac)

Compute the Jacobian. jac[i] is d`mdata`/d`pars[i]`.

extract(pars, perr, cov)

Extract fit results into the object for ease of inspection.

finalize_setup()

If the component has subcomponents, this should set their name, setguess, setvalue, and setlimit properties. It should also set npar (on self) to the final value.

model(pars, mdata)

Modify mdata based on pars.

prep_params()

This should make any necessary calls to setvalue or setlimit, though in straightforward cases it should just be up to the user to do this. If the component has subcomponents, their prep_params functions should be called.