1. Parameter Estimation
PyDSTool provides a thin wrapper around the SciPy versions of the Levenberg-Marquardt least squares optimizer and the one-parameter BoundMin bounded minimizer, under a unifying class type ParamEst, which is implemented in PyDSTool/Toolbox/ParamEst.py. This class provides basic features such as a default residual function that is compatible with PyDSTool Model objects, easing the changing of model parameter values, computation of trajectories, and setting of additional parameters for each call to the function.
1.1. Basic usage
The user begins by creating and initializing a ParamEst class instance with a Model system, the number of free parameters in the model to be fitted, the dependent variables names to use, and the residual function to use. Either integration initial conditions or model parameters can be used as the free parameters.
Residual functions must have a signature of the form (s, p, tmesh [, <optional additional argument sequence>]). The first argument s refers to the ParamEst class instance initialized by the user, which provides the function access to the model and other user-specified settings. The p argument will be an array of new parameter values, given in the same order as s.freeParNames.
An explicit Jacobian defined in a model system (such as that calculated using the Symbolic toolbox differentiation functions) will automatically be used by the class.
Once the ParamEst object has been created, its run method can be called in order to execute the actual optimization algorithm. Arguments to this call allow additional arguments to be sent to the residual function, and algorithmic tolerances to be set. By default, the ParamEst class calls the optimizers with their 'full output' options set, and stores the results in the dictionary s.pestResult, as well as returning a copy of that dictionary. Note that this dictionary will be overwritten on the next call to run.
One benefit of organizing a parameter estimation problem inside a class such as ParamEst is that it gives the residual function greater power to store state inside its associated class instance.
Example of a minimal setup for features and model context from /tests/joe_pest.py (version 0.87):
# ---- make features for measuring similarity
# feature residual tolerance
ftol=3e-3
# use L2-norm of data (sum of squares)
L2_similarity = L2_feature('L2_similar', pars=args(t_samples=tmesh,
trange=tdomain,
coord='v',
tol=ftol))
# condition for feature is that it is present in data (True)
# (could also specify that it is not present using False)
pest_condition = condition({L2_similarity : True})
# trivial sub-classes of external and internal model interfaces
# that will interact with each other
class ext_iface(extModelInterface):
# holds the data (external from the model)
pass
class int_iface(intModelInterface):
# holds the test model
pass
# embed the reference trajectory and measurement condition/features
# into an instance of the external interface
pest_data_interface = ext_iface(reftraj,
pest_condition)
# generate a context for the model between the external interface
# instance and the internal interface class (not an instance yet)
c = context([ (pest_data_interface, int_iface) ])
# specify which parameters we will fit
est_parnames = ['gl', 'vl']
# parameter estimation
print 'Starting Least Squares parameter estimation'
print 'Goal pars are gl = ', est_pars_ref['gl'], ' vl = ', est_pars_ref['vl']
pest_pars = LMpest(freeParams=est_parnames,
testModel=testModel,
context=c
)
pestData_par = pest_pars.run(parDict={'ftol':ftol,
'xtol':1e-3},
verbose=True
)
For further information see examples in PyDSTool/tests such as pest_test1.py for initialization of ParamEst class instances, overriding of default residual functions, and running of the optimizers.
1.2. Algorithms available
1.2.1. Levenberg-Marquardt least squares (class LMpest)
This works in an N-dimensional parameter space on an n-dimensional model trajectory. This sub-class wraps the leastsq function from scipy.minpack. This algorithm is inherently unconstrained, but can be adapted to incorporate soft constraints in the residual function (see the provided example scripts).
The class attribute _memory, which is pre-declared and provided as a memory for the residual functions between calls. It is declared to be the array zeros(self._tmesh_len*self._depvar_len, dtype=float64), to provide temporary storage while building the return array of residuals (see residual_fn_LMpest_default_Ndim for an example).
1.2.2. Bounded minimization (class BoundMin)
This works only in a one-dimensional parameter space for a one-dimensional model trajectory. This sub-class wraps fminbound from scipy.optimize. It minimizes subject to bound constraints.
1.3. Qualitative fitting
A new feature to v0.85 is to be able to fit the extrema of a trajectory with those of data using conventional optimization algorithms, such as least squares. This qualitative fitting approach is under development but appears to work well in instances where standard approaches for which the fitting problem is ill-conditioned by virtue of having a "fitness landscape" that is too flat near typical initial conditions.
A robust way to fit extrema is to define events in the test model to find them accurately. The tools provided here assume that you have added appropriate events to your model to do this. This is described further in the following docstring for the example residual function residual_fn_LMpest_extrema(s, p, tmesh) exported from ParamEst.py:
Example residual function for the qualitative fitting of local
extrema locations in model trajectory to those of data, using least squares. This assumes that dictionaries containing the extrema data have been pre-computed and placed in alg_params, as under keys 'data_mins_t', 'data_mins_v', 'data_maxs_t', 'data_maxs_v', for the time and variable value locations of the maxima and minima. Each dictionary of extrema data is expected to be an array of values keyed by coordinate name.
The generator specified by gen_name should exist in the model s.testModel given to the ParamEst object at creation. This generator should contain separate events for detecting maxima and minima in each coordinate being fitted (e.g. zero crossing of first derivative of a coordinate, in appropriate direction). For a coordinate named 'x' these events are expected to have the names 'min_ev_x' and 'max_ev_x'.
num_expected_mins and num_expected_maxs is per variable.