iFit: iData objects fitting


  1. Fitting a model to the Data
    1. Getting most of the fit results (4th output argument)
  2. Specifying/configuring the optimization method
    1. Choosing an optimization method
    2. Configuring the optimization method (4th input argument)
    3. Monitoring the performance of the optimization method and the fit process
  3. Model parameter constraints
    1. Fixed parameters
    2. Parameters varying within limits
    3. Limiting parameter change
  4. Estimating the model parameter uncertainties and the fit quality
  5. Building model functions
  6. 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.

fits(iData)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:
  1. a: an iData object, or an array of objects
  2. model: the name of the function or a function handle
  3. 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:
  1. parameters: is the best fit parameter set for the model obtained during the optimization
  2. criteria: is the criteria value (least square χ² by default - see below on how to change this setting)
  3. message: is the final optimization routine state (converged, failed, ...)
  4. 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:
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:

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 ])
fits with options.OutputFcn=fminplot iData fits(gauss) with plot
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:
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

ifitmakefuncThe 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 ILL, Grenoble, France <www.ill.eu>