EasySpin: esfit

esfit

Synopsis

Least-squares fitting of EPR spectra.

esfit(SimFnc,expSpc,Sys0,Vary,Exp);
esfit(SimFnc,expSpc,Sys0,Vary,Exp,SimOpt);
esfit(SimFnc,expSpc,Sys0,Vary,Exp,SimOpt,FitOpt);
BestSys = esfit(...)
[BestSys,BestSpc] = esfit(...)

See also the tutorial on fitting.

Description

esfit fits EPR spectra simulated with garlic, pepper, and chili to experimental spectral data using least-squares fitting algorithms.

The input parameters are the following:

SimFnc is the name of the simulation function, either 'garlic', 'pepper', 'chili', 'saffron' or 'salt'.
expSpc is a 1D array containing the experimental spectral data.
Sys0 is a spin system structure containing all the parameters of the spin system. For parameters that are to varied during the fitting process, the values in the structure represent the centers of the search ranges.
Vary contains the approximate maximum variations for the spin system parameters allowed in the fitting. If a spin system parameter should be included in the fitting, give it a non-zero value in this structure. If it should not be included, set its value in this structure to zero, or don't include it at all.
```
Sys0.lwpp = 5;
Vary.lwpp = 2;    % the linewidth is searched over a range of about 3 to 7
Vary.lwpp = 0;    % the linewidth kept constant at the value given in Sys0
```
SimOpt contains simulation options that are forwarded directly to the simulation function. Refer to the documentation of chili, garlic, pepper, and salt for details.
FitOpt is a structure containing settings for the optimization algorithms used in esfit. The possible settings are discussed further down on this page.

If no output is requested, esfit plots the result of the fitting. Several optional outputs are available:

BestSys contains the spin system parameters of the best fit.
BestSpc contains the simulated spectrum that gives the best fit.

The structure FitOpt contains fitting options. The most important ones select the fitting algorithm and the function to use for computing the residuals:

Method
A string indicating the least-squares fitting method to be used, consisting of up to two keywords. One keyword specifies the algorithm, and another one selects the way residuals are computed. Some examples:
FitOpt.Method = 'simplex fcn'; FitOpt.Method = 'genetic int';

The keywords for choosing one of the available algorithm are

'simplex': Nelder-Mead simplex algorithm
'levmar': Levenberg/Marquardt
'montecarlo': Monte Carlo random search
'genetic': Genetic algorithm
'grid': Systematic grid search
If you don't specify anything, the default 'simplex' is used.
The second keyword specifies what the residuals and consequently the root-mean-square deviation (rmsd) should be computed from. Usually the residuals are computed by taking the difference of the experimental and the simulated spectrum. However, the spectra can be transformed before computing the residuals.

'fcn' means take the experimental and simulated spectral data directly.
'int': integrate the two spectra.
'dint': use the double integrals of the two spectra.
'diff': compute the derivatives of the two spectra.
'fft': FFT both spectra before computing the residuals.
The default value is 'fcn'. If there are many resolved lines, 'int' gives better convergence.

Scaling
This string specifies the scaling method. There are different ways a simulated spectrum can be scaled and shifted to overlap as much as possible with the experimental one. The possible values are

'maxabs': make the maximum absolute of simulated and experimental spectrum equal
'minmax': rescales and shifts so that minimum and maximum of the two spectra are identical
'lsq': obtains the scaling factor by least-squares fitting, without shifting
'lsq0': same as 'lsq', but includes a constant baseline correction
'lsq1': same as 'lsq', but includes a linear baseline correction
'lsq2': same as 'lsq', but includes a quadratic baseline correction

If the baseline of the experimental spectra has been corrected beforehand, the best choice is 'maxabs'. This is also the default value. If the noise level is high, 'lsq' is a better choice:
FitOpt.Scaling = 'lsq';
To include corrections for a baseline in the expermental spectrum, use 'lsq0, 'lsq1, or 'lsq2.

There are some more options which are recognized by all fitting algorithms:

maxTime Time, in minutes, after which esfit will terminate even if the fitting has not yet converged.

PrintLevel A number, 0 or 1. If set to 0, the fitting functions don't write anything to the command window. If set to 1 (default), they print information about the progress of the fitting.

Plot If set to 1, the fitting progress is plotted. If set to 0, all plotting is suppressed.

Each of the available fitting algorithms can be fine-tuned using a set of parameters.

Parameters for Nelder/Mead downhill simplex:

RandomStart If set to 1, the starting point in the parameter space is chosen randomly, within the given limits. If set to 0, the center of the parameter space is used as starting point. By default, RandomStart is 0.

TolStep Termination tolerance for the length of the parameter step. This number refers to the rescaled fitting parameters, as described above.

TolFun Termination tolerance for function value change.

SimplexPars An array of four elements [rho chi psi sigma], where rho is the reflection coefficient, chi is the expansion coefficient, psi is the contraction coefficient, and sigma is the shrink coefficient. The default values are [1,2,0.5,0.5].

delta Size parameter of the initial simplex. The default value is 0.1.

Parameters for Levenberg-Marquardt:

RandomStart If set to 1, the starting point in the parameter space is chosen randomly, within the given limits. If set to 0, the center of the parameter space is used as starting point. By default, RandomStart is 0.

TolStep Termination tolerance for the length of the parameter step. This number refers to the rescaled fitting parameters, as described above.

TolFun Termination tolerance for function value change.

lambda Starting value of Marquardt parameter λ, default value is 0.001.

delta Step size for computing the finite-difference approximation of the Jacobian. Default is 1e-7.

Parameters for Monte Carlo:

N Number of random trial simulations, default 200.

TolFun Termination tolerance for error function value change.

Parameters for the genetic algorithm:

PopulationSize A number giving the size of the population, that is the number of parameter sets and simulations in one generation. The default value is 20, but for fittings with many parameters, this value should be increased.

maxGenerations A number specifying the maximum number of generations the algorithm should run. After this number has been reached, the algorithm terminates, no matter how good or bad the best fit so far is. The default value is 40, but has to increased for fittings with many parameters. If only very few parameters are fitted, maxGenerations can be decreased.

Parameters for the grid search:

GridSize A number or an array that specifies how many grid points there should be for each parameter. If one number is given, it is valid for all parameters. For example, let's assume that one g value and the linewidth are being fitted:
Vary.g = [0 0.001 0]; Vary.lw = 0.2;
Then GridSize can contain 1 or 2 numbers:
FitOpt.GridSize = 10; % 10 points for each parameter, making 100 grid points total FitOpt.GridSize = [20 3]; % 10 points along g and 3 along lw, giving a total of 60

Algorithms

Internally, esfit rescales the parameters to be fitted. E.g. if the linewidth should be fitted and is expected to be around 5 mT, plus minus about 1.5 mT, Sys.lw = 3; and Vary.lw = 1.5. After rescaling, esfit treats 5 mT as x = 0, 3.5 mT as -1 and 6.5 mT as +1.

For details about the algorithms implemented in esfit, see e.g.

Press et al, Numerical Recipes: The Art of Scientific Computing, 3rd edition, Cambridge University Press, 2007, Chapters 10 and 15.
Spalek et al, J.Chem.Inf.Model 2005, 45, 18-29.

Examples

Here is a very simple example. Let's assume the experimental data are stored in expspec. Then the following code performs a least-squares fitting using the simplex algorithm.

Exp.mwFreq = 9.5;
Sys0.g = [2.1 2.2]
Sys0.lwpp = 0.1;
Vary.g = [0.05 0.02];
esfit('pepper',expspc,Sys0,Vary,Exp);

See the example section for a full example. Also, read the tutorial about fitting.

`maxTime`	Time, in minutes, after which `esfit` will terminate even if the fitting has not yet converged.
`PrintLevel`	A number, 0 or 1. If set to 0, the fitting functions don't write anything to the command window. If set to 1 (default), they print information about the progress of the fitting.
`Plot`	If set to 1, the fitting progress is plotted. If set to 0, all plotting is suppressed.

`RandomStart`	If set to 1, the starting point in the parameter space is chosen randomly, within the given limits. If set to 0, the center of the parameter space is used as starting point. By default, `RandomStart` is 0.
`TolStep`	Termination tolerance for the length of the parameter step. This number refers to the rescaled fitting parameters, as described above.
`TolFun`	Termination tolerance for function value change.
`SimplexPars`	An array of four elements `[rho chi psi sigma]`, where `rho` is the reflection coefficient, `chi` is the expansion coefficient, `psi` is the contraction coefficient, and `sigma` is the shrink coefficient. The default values are `[1,2,0.5,0.5]`.
`delta`	Size parameter of the initial simplex. The default value is 0.1.

`N`	Number of random trial simulations, default 200.
`TolFun`	Termination tolerance for error function value change.

`PopulationSize`	A number giving the size of the population, that is the number of parameter sets and simulations in one generation. The default value is 20, but for fittings with many parameters, this value should be increased.
`maxGenerations`	A number specifying the maximum number of generations the algorithm should run. After this number has been reached, the algorithm terminates, no matter how good or bad the best fit so far is. The default value is 40, but has to increased for fittings with many parameters. If only very few parameters are fitted, `maxGenerations` can be decreased.