iFit: iData objects fitting
- Fitting
a
model
to
the
Data
- Getting most of the
fit results (4th output argument)
- Specifying/configuring the optimization
method
- Choosing an
optimization method
- Configuring the
optimization method (4th input argument)
- Monitoring the
performance of the optimization method and the fit process
- Model parameter constraints
- Fixed parameters
- Parameters varying
within limits
- Limiting parameter
change
- Estimating the model parameter
uncertainties and the fit quality
- Building model functions
- Specifying the optimization criteria
Commands we use in this page: iData/
fits,
ieval iFuncs models and
iOptim
routines
This documentation details the procedure that can be used to find
optimal parameter set from a model in order to match (fit) an iData
object Signal with axes.
fit: { vary vector
p ∈
ℜd
so that
distance(Signal, Model) is
minimal }
where
distance is the criteria
for the fit, which measures how far the model is from the data. There
is a number of possible distance definitions, detailed
below.
Fitting
a
model
to
the
Data
The fits
function for iData objects provides a simple
mean to find the best model parameters which fits a data set.
>> a=load(iData, [ ifitpath 'Data/sv1850.scn' ])
>> p=fits(a);
'Amplitude' 'Centre' 'HalfWidth' 'Background'
0.6786 1.0008 0.0035 0.0002
will fit the data Signal/Monitor
to a default 1D Gaussian function 'gauss',
starting
from
automatically
guessed
initial
model
parameters.
Of
course,
it
is
possible
to
chose explicitly
the model function, and the starting parameter set with e.g.
>> p=fits(a, 'gauss'); % specify model function to use: gauss, and use guessed parameters
>> p=fits(a, @gauss, [ 0.5 1 0.01 0 ]); % specify the starting parameters for the model function
The result of the fits method
is the model parameter set that matches best the data. The number of
fit parameters is not limited.
A number of model functions is available from the iFuncs
sub-library. The list of all
available fit functions can be obtained from the command:
>> fits(iData);
Refer
to the section below in order to learn
how to define new fit models with the ifitmakefunc
tool. You can get information about a model by
inquiring its name. When a single input parameter is specified,
starting parameters
are guessed :
>> gauss
>> gauss('identify') % return informations about the function
>> gauss('plot') % plot the function with a default parameter set
>> gauss(p, x) % evaluate the function with parameters 'p' and axis 'x'
>> gauss(p)
ans =
Type: 'iFit fitting function'
Name: 'Gaussian (1D)'
Parameters: {'Amplitude' 'Centre' 'HalfWidth' 'Background'}
Dimension: 1
Guess: [1 6.2172e-17 0.5861 0.1000]
Values: [1x100 double]
Axis: {[1x100 double]}
A fast way to define a new model and fit it to the data is to directly
execute, with p the parameter vector, and x,y,... the axes:
>> p=fits(a, 'p(1)*x+p(2)');
which will transparently call the ifitmakefunc
tool to create the model function, and then use it for the fit.
In
order
to
estimate
the
values
of
the
model
for
the
fit
parameters
onto
the iData object data set axes, one can use the ieval method,
which is the same as Matlab feval, but for iData
objects:
>> b= ieval(a, 'gauss', p) % evaluate the 'gauss' model onto the 'a' axes
>> plot([ a b ]) % with parameters 'p', and plot both data and fit
The evaluated model contains Parameter
and Model aliases that hold
the model parameters and description.
The fits and ieval functions work with any data
and model dimensionality, and can also build multi-dimensional models
from lower dimensionality functions (e.g. making a 2D Gaussian from
multiplied perpendicular 1D Gaussians).
Getting most of the
fit results (4th output argument)
The fits method can return
additional arguments
>> [parameters,criteria,message,output]= fits(a, model, initial_parameters,...)
The input arguments allow to specify:
- a: an iData object, or
an array of objects
- model: the name of the
function or a function handle
- initial_parameters: a
set of starting parameters for the model. Using an empty set will trigger an automatic estimate (guess) of the
starting model parameters.
The returned arguments are similar to those returned by the
Matlab fminsearch non-linear
optimizer:
- parameters: is the best
fit parameter set for the model obtained during the optimization
- criteria: is the
criteria value (least square
χ² by default - see below on how to
change this setting)
- message: is the final
optimization routine state (converged, failed, ...)
- output: is a structure
that holds all the information gathered during the fit procedure
The most interesting output argument result is the 4th one output,
which
provides significantly more information such as:
- output.parsBest: best
fit parameter set for the model
- output.criteriaBest: best criteria value for the model
- output.modelValue: Last model
evaluation (iData), that is ieval(a, model, parameters)
- output.algorithm: Algorithm/solver
used
(char)
- output.message: Message
which details the final state of the optimizer (char)
- output.funcCount: Number of
function evaluations performed during the fit (double)
- output.iterations: Number of
iterations performed during the fit (double)
- output.parsHistory: Parameter
set history during optimization (double array)
- output.criteriaHistory: Criteria
history during optimization (double vector/array)
- output.parsHistoryUncertainty: uncertainty
on fit parameters (half-width) obtained from the optimization history
(vector)
- output.parsHessianUncertainty: uncertainty
on fit parameters (half-width) obtained from the Hessian matrix, which
is extremely sensitive to noise. Prefer the output.parsHistoryUncertainty value
(vector)
- output.corrcoef: fit
correlation coefficient (closer to 1 indicates a good fit), as obtained from Matlab corrcoef.
- output.modelName: model
name
- output.modelInfo: model
information (structure)
- output.parsNames:
parameter names from the model
Last, the starting parameter 'guess' can also be specified as a
structure which fields contain numbers, or as a string with members separated with the ';' character, such as in the
following example:
>> p.Amplitude=0.5; p.Center=1; p.HalfWidth=0.0035; p.Background=1e-4; % optimize named parameters
>> fits(a,'gauss', p) % The result is also returned as a structure.
>> fits(a,'gauss','Amplitude=0.5; Center=1; HalfWidth=0.0035; Background=1e-4') % create the structure above and fit...
All model named parameters must be specified.
Specifying/configuring the optimization
method
The iData fits
method is a wrapper to any of the iOptim
optimizers,
with a specified fit criteria. Each optimizer can be
customized (e.g. the maximum
number of iterations, stop/convergence conditions, ...) with an options structure.
>> [parameters,criteria,message,output]= fits(a, model, initial_parameters, options)
Choosing an
optimization method
The iOptim sub-library provides 21 different
optimization techniques. No doubt that choosing one at first sight is
difficult. We provide below a set of preferred optimizers, based on a
careful comparison explained in the iOptim
documentation.
Optimizer
|
Description
|
Mean
Success
ratio
(%)
|
fminpso |
Particle Swarm
Optimization |
97.9 / 84.1 |
fminimfil |
Unconstrained Implicit filtering
|
93.5 / 53.8
|
fminralg |
Shor R-algorithm (avoid when noisy) |
88.9 / 16.2
|
fminhooke |
Hooke-Jeeves direct search |
97.1/ 56.3
|
fmincmaes |
Evolution Strategy with
Covariance Matrix Adaptation |
88.9 / 71.2
|
fminsimpsa |
simplex/simulated
annealing |
97.1 / 84.8
|
fminsce |
shuffled complex evolution
|
95.7 / 85.2 |
fminpowell |
Powell with Coggins line search
|
99.2 / 51.7 |
Table
1:
A
selection
among
the
most
efficient
and
fast
optimization
methods. The two success ratio are given for continuous and noisy functions resp.
The choice of the optimizer is done through a 4th input argument options
to fits (see below).
The list of all available optimizer
methods can be obtained from the command:
>> fits(iData)
Configuring the
optimization method (4th input argument)
In order to use these optimizers, one just has to specify their names
as
a fits input argument, e.g. options='fminimfil' :
>> [parameters,criteria,message,output]= fits(a, model, initial_parameters, 'fminimfil')
will select the Unconstrained
Implicit filtering coupled with BFGS fminimfil
optimizer to perform the fit. The default optimizer configuration will
be used, as obtained from the optimset
function, or from the method itself with 'defaults' as parameter. This
latter call (i.e. not with optimset) provides more informations:
>> options=fminimfil('defaults') % get the default optimizer configuration parameters
The optimizer configuration is
a structure which members enable to tune the behavior of the
optimization process. Each of the fields can be changed, e.g.
>> a=load(iData, [ ifitpath 'Data/sv1850.scn' ])
>> options=fminimfil('defaults')
>> options.TolFun=0.01;
>> p=fits(a, model, [], options); % fit with the customized options, and guessed starting parameters
In addition to the default optimset
structure fields, the iData fits
and the iOptim optimizers use
additional members in the structure:
- options.Display: Level
of display [ off | iter | notify | final ]. Default is 'off'
- options.MaxFunEvals: Maximum number of
function evaluations allowed
- options.MaxIter: Maximum number of
iterations allowed
- options.TolFun: Termination tolerance on
the function value (absolute value or change). Can be specified in '%x' for a relative value.
- options.TolX: Termination tolerance on
parameter change. Can be specified in '%x' for a relative value.
- options.OutputFcn: Name of an output
function (see below). You may use 'fminplot', which is
provided in iOptim.
This function should have the syntax fminplot(pars,
optimValues, state).
- options.FunValCheck: Check for invalid
function/model values, such as NaN or complex and abort.
- options.algorithm: is a description of the
optimization method (solver)
- options.optimizer: is the name of the
optimization function (solver)
- options.criteria: the
criteria to use with syntax criteria(Signal,
Error, Model). The default is least-square error (see below).
- and more, depending on the optimizer method and implementation
Monitoring the
performance of the optimization method and the fit process
In order to follow the optimization process, you may define a call to a
user function at each optimization iteration. A default plotting
facility has been implemented as the fminplot
function:
>> a=load(iData, [ ifitpath 'Data/sv1850.scn' ])
>> options=fminimfil('defaults')
>> options.OutputFcn='fminplot';
>> p=fits(a, 'gauss', [], options)
p =
0.6263 1.0008 -0.0037 0.0002
>> b = ieval(a, 'gauss', p)
>> figure; plot([ a b ])
Left: The
options.OutputFcn='fminplot' window
on
the
right,
showing
the
criteria
evolution
with
the
optimization
iteration
and
up
to the 3 first fit parameters. The
red
dot indicates the
current/final parameter set and criteria.
Right: The
final fit result on the left.
Also, as explained
above, it is possible
to obtain, on optimization completion, the whole criteria, and
parameter set history, as a function of iterations/criteria evaluations.
>> [p,criteria,message,output]= fits(a, 'gauss', [], options)
>> output.parsHistory
>> output.criteriaHistory
Model parameter constraints
In many cases, the model parameters are to be constrained. Some optimization specialists call these restraints, that is parameter values constraints. This
includes some fixed values, or bounds within which parameters should be
restricted. This is given to the fits
method by mean of a 5th input
argument constraints
:
>> [parameters,criteria,message,output]= fits(a, model, initial_parameters, options, constraints)
Fixed parameters
To fix some of the model parameters to their starting value, you just
need to define constraints as
a vector with 0 for free parameters, and 1 for fixed parameters, e.g. :
>> a=load(iData, [ ifitpath 'Data/sv1850.scn' ])
>> p=fits(a, 'gauss', [], 'fminimfil', [ 1 0 0 0 ])
** Minimization performed on parameters:
'Amplitude' 'Centre' 'HalfWidth' 'Background'
0.5936 0.9998 0.0018 0.0050
p =
0.5936 1.0008 -0.0037 0.0002
will fix the first model parameter, which is here the Amplitude. This parameter name
can be checked
by simply entering the name of the model, which returns some
information structure :
>> gauss
ans =
Type: 'iFit fitting function'
Name: 'Gaussian (1D)'
Parameters: {'Amplitude' 'Centre' 'HalfWidth' 'Background'}
Dimension: 1
Guess: [1 6.2172e-17 0.5861 0.1000]
Values: [1x100 double]
Axis: {[1x100 double]}
A similar behavior is obtained when setting constraints as a structure
with a fixed member :
>> constraints.fixed = [ 1 0 0 0 ];
>> p=fits(a, 'gauss', [], 'fminimfil', constraints);
The constraints.fixed vector should
have the same length as the model parameter vector.
Parameters varying
within limits
If one needs to restrict the exploration range of parameters, it is
possible to define the lower and upper bounds of the model parameters.
This can be done by setting the 5th fits
argument to the lower bounds lb,
and the 6th to the upper ub,
e.g. :
>> a=load(iData, [ ifitpath 'Data/sv1850.scn' ])
>> p=fits(a, 'gauss', [], 'fminimfil', [ 0.5 0.8 0 0 ], [ 1 1.2 1 1 ])
lb ub
A similar behavior is obtained by setting constraints as a structure
with members min and max :
>> constraints.min = [ 0.5 0.8 0 0 ];
>> constraints.max = [ 1 1.2 1 1 ];
>> p=fits(a, 'gauss', [], 'fminimfil', constraints);
The constraints min and max vectors
should have the same length as the model parameter vector.
Limiting parameter
change
Last, it is possible to restrict the change rate of parameters by
assigning the constraints.steps
field to a vector. Each non-zero value then specifies the absolute
change that the corresponding parameter can vary between two optimizer
iterations.
In short, the constraints structure can have the following members,
which all should have the same length as the model parameter vector:
- constraints.fixed: 1
for fixed parameter, 0 for free parameters
- constraints.min: the
minimum value for each parameter. -Inf is supported.
- constraints.max: the
maximum value for each parameter. +Inf is supported.
- constraints.steps: the
maximum change between iterations for each parameter. +Inf is supported.
All these constraints may be used simultaneously. NaN values in these constraints are ignored (the corresponding parameters are not constrained).
Estimating the model parameter uncertainties
and the fit quality
Some theoretical notes about the goodness of fit are indicated in
the iOptim help. The
uncertainty on the fit parameters can be obtained from output.parsHistoryUncertainty and output.parsHessianUncertainty. Additionally,
the output.parsHessianCorrelation
indicates the cross-correlations between parameters (off-diagonal
values larger that 0.7 indicate cross-correlated parameters).
The fit quality can be assessed from the output.corrcoef returned value
which goes from 0 (very bad fit) to 1 (perfect fit). A value higher
than 0.95 is usually good.
Building model functions
The
ifitmakefunc
tool has been designed to automatically create model functions from
either a single expression, or a more detailed description. Refer to
the iFuncs/Model Builder.
>> h=ifitmakefunc; % pops-up a dialog to define the new fit function/model
>> h=ifitmakefunc('p(1)*exp( (x-p(2))/p(3) )');
>> fits(a, h)
The resulting function has the ability to identify itself ('identify', provide
detailed informations), compute automatic starting parameters ('guess'), display
itself ('plot'), and
evaluate its value of course. It can be directly used with
ieval and fits, either with their name, or
their function handle.
It is even possible to directly call the fitting method and create the
model function on the fly, which makes the fit much easier for simple
functions that can be written as a single expression:
>> fits(a, 'p(1)*exp( (x-p(2))/p(3) )'); % does the same as above, in a single command
To assemble existing functions into new ones, you may use e.g.:
>> h=ifitmakefunc('gauss(p(1:4),x)+lorz(p(5:8),x)','p(8)=0;');
which creates a new function which is the sum of a Gaussian and a
Lorentzian. The second redundant Lorentzian Background p(8) parameter
is forced to 0 so that it does not correlate with the Gaussian
Background p(4). Other
function information/parameter names (here not specified) are
guessed/defaulted.
It is even possible to convolute and correlate functions and data sets
as new function definitions. Refer to the iFuncs page.
In case the model requires additional arguments, just pass them to the fits method (arguments above the 5th)
>> p=fits(a, model, [ 0.5 1 0.01 0 ],'','',additional_arguments);
assumes the model function has syntax
model(p, axes from object, additional_arguments)
Specifying the optimization criteria
The default optimization criteria is the 'least square error':
χ² = ∑ (Signal-Model)².
When the Monitor is defined (see iData
definition), the Signal is normalized, so that
χ² = ∑
(Signal/Monitor-Model)².
When the Error on the Signal is available, the weighted criteria reads
χ² = ∑
(Signal/Monitor-Model)²/Error².
The options.criteria can be defined as the name of a criteria function
with syntax criteria(Signal, Error,
Model) where all arguments should have the same dimension, and
the default criteria is:
options.criteria
=
'least_square';
% internal 'fits' least square definition
The Signal will be normalized to the Monitor prior to calling the
criteria. The Error is used as weight.
The criteria is then normalized to the number of degrees of freedom by dividing it by
χ² ← χ² /(number_of_data_points_in_the_Signal - number_of_parameters -1)
The following pre-defined functions can be used as criteria:
criteria(Signal, Error,
Model) |
Expression
|
comment
|
least_square
|
(|Signal-Model|/Error)2
|
non-robust. Ref.
|
least_absolute
|
|Signal-Model|/Error
|
robust. Ref.
|
least_median
|
median(|Signal-Model|/Error)
|
robust, scalar. Ref.
|
least_max
|
max(|Signal-Model|/Error)
|
non-robust, scalar. Ref.
|
E.
Farhi - iFit/iData fitting -
$Date: 2012-02-11 02:27:24 $ $Revision: 1.24 $
- back to
Main iFit Page