FitError

FitError

The FitError function is used for estimating error of the fitting parameters. Typically, it should be run using the same fitting settings as the initial fit was acquired with, in order to obtain accurate error analysis. The function takes a spectrum fit and the calculated spectrum from that fit, adds noise to it, and refits it. This is done many times (~100s), in order to estimate the effect that experimental noise has on fitted parameters. Statistics of the refit parameters can be used to estimate statistics of the experimental parameters. Note that application of this method is limited, however. First, it assumed that noise in the time domain is uncorrelated, white noise, so that error induced by spectrum artifacts, incomplete phase cycling, etc. will not be accounted for. Second, the error values produced are the representative of experimental error only if the fit is a good representation of the ‘correct’ values of the peak parameters. Therefore, in regions where there are missing peaks, or if a peak is picked which is actually noise, then the error in those situations will likely not represent the real error (one can consider the ‘FitError’ function to be simulating experimental repetition). This is discussed in more detail in the accompanying publication.

The FitError function does not repeatedly refit the full spectrum. Rather, for each peak that is evaluated for error, FitError fits a partial spectrum around that peak, using the PartialFit function. FitError requires that if the error of a peak is determined from a partial fit, then that partial fit must also include neighboring peaks that fall within a certain range of the peak of interest (see ‘par.NB_cutoff’ below). Also, the region fitted around the peak and its neighbors must have a certain width (see ‘par.H_cutoff’). However, since now some sections must contain multiple peaks, it is additionally possible to analyze the error of those peaks- if the section can be modified so that these additional peaks also have their neighbors and a wide enough range around the peak. This set of peaks is refererred to as a group in INFOS. To prevent a group from becoming too large (and therefore slow to fit), FitError restricts the number of peaks and the frequency range (‘par.max_peaks’ and ‘par.max_width’). These restrictions mean that the fitting of individual groups can be performed without breaking the spectrum into further sub-spectrum, eliminating the need for multiple iterations of fitting. Because the number of parameters is limited, the fitting is relatively fast. Note that it is then possible that a peak will be fitted in a group, but it is not possible to include its neighbors, so that error analysis is not performed on that peak. In this case, that peak will have its own group, in order to obtain statistics. This will lead to redundancy in fitting some peaks, but in spite of this, breaking the spectrum into groups is significantly faster than repeating the fit of the complete spectrum many times. After fitting a spectrum, FitError can be called as follows

fit=FitSpec(spec,fitpar);
err=FitError(fit,par);

It is also possible to call without any parameters (one or both of the fitting parameter arguments can be left out.

fit=FitSpec(spec);
err=FitError(fit);

The ‘fit’ structure is simply the output of the FitSpec function. The ‘par’ structure can be constructed as for FitSpec. However, note that the values in ‘par’ used during the initial fit are passed to the FitError function via the substructure ‘fit.par’. Since one should use the same fitting parameters in the initial fit and in the FitError analysis, it typically is not a good idea to pass new parameters via the ‘par’ structure. The exception is that a few additional parameters are used by FitError for determining number of refits and for forming the groups of peaks to be fit. Note these parameters all have default values, and so do not need to be provided.

par.N: This is the number of times the refitting should be performed. Default is 100, which gives reasonable estimates of the standard deviation in most cases. However, if a full distribution of a parameters is required (for example, if distribution is not a normal distribution), then many more (~1000) iterations may be needed.

par.NB_cutoff: This field is determines how far to search around a peak for its neighbors. For a peak of interest, ‘par.NB_cutoff’ is multiplied by that peak’s height. Then, in each dimension, FitError determines what distance from the peak center the peak falls below this amplitude. Neighboring peaks are then defined as any peak falling within this range in all dimensions. Default is 10% of the maximum peak height (‘par.NB_cutoff=0.1’).

par.H_cutoff: This field determines how wide the fitted area around peaks in a group needs to be. The value is defined similarly to ‘par.NB_cutoff’, with the width required being determined from a fraction of the peak height. Then, a group must be wide enough that all peaks in the group have at least the width around them defined by ‘par.H_cutoff’.

par.max_width This is the maximum spectral width a group may have. This must have one entry per spectrum dimension. Note that when a group is formed, it is possible that the neighbors around the first peak in the group are spread out enough to already exceed the maximum width. In this case, ‘par.max_width’ will be overridden, although no further peaks will be added to such a group. Default is 1.5x the size specified by ‘par.grid’ in the initial spectrum fit.

par.max_peaks: This is the maximum number of peaks allowed in a group. As with ‘par.max_width’, if the first peak has more neighbors than allowed by ‘par.max_peaks’, then this parameter will be overridden. Default is 10.

par.index: This allows the user to determine which peaks should have error analysis performed. Field is the same that for PartialFit (section 3.6.3). It should either be a logical index with the same size as ‘fit.I’, or a vector listing which peaks should be analyzed for error.

The output of FitError then contains statistics calculated for the individual fitting parameters, in addition to a record of the fitting results of all ‘par.N’ fitting iterations. These are stored in the following fields

err.par: Contains values used for error analysis as described above (‘par.N’, ‘par.NB_cutoff’, etc.)

err.fit: Returns the initial fit (simply the input variable returned)

err.I_std: Standard deviations of the peak amplitudes

err.deltaX_std Standard deviations of the peak position in the Xth dimension

err.lwX_std Standard deviations of the peak linewidth in the Xth dimension

err.FWHMX_std Standard deviations of the peak FWHM in the Xth dimension

err.int_std Standard deviations of the peak integral in the Xth dimension

err.groups Cell with the results of fitting each of the groups (has same number of elements as the number of groups). Inside each cell of ‘err.groups’, the results of the ‘par.N’ refittings are stored. Therefore, this contains all the same fields as in the usual ‘fit’ structure (section 2.4), excepting that the spectra (‘fit.spec’, ‘fit.resid’) have been removed to reduce the variable size. However, each fit parameter (‘fit.I’, ‘fit.delta1’, etc.) is now an Nxpar.N matrix, with N being the number of peaks in the group, and par.N the number of repetitions. There is also an additional field, ‘err.groups{k}.calc_list’ which lists for each peak in the group its index in the original fit. Note that the field ‘err.groups{k}.par’ is as described in section 2.2, but also includes fields ‘err.groups{k}.par.fit_range’ which indicates the range of the spectrum that has been fitted, and ‘err.groups{k}.par.index’ which indicates which peaks are fitted in a group.

err.group_index This list indicates which group contains each peak. Note that some peaks can be contained in more than one group, but statistics should be analyzed from the group indicated here.

err.subgroup_index This list indicates for each peak, which element of the group corresponds to this peak

To understand indexing of groups better, a quick example on how to determine the histogram of a peak’s amplitude is shown. Suppose that ‘x’ is an edges vector for a histogram (see histc.m), and we want the histogram of the peak with index ‘k’. Then, the histogram is given by

 h = histc(err.groups{err.group_index(k)}.I(err.group_subindex(k),:),x-diff(x(1:2))/2);

One sees that the correct group was found with the index ‘err.group_index(k)’, and the correct element of that group was found with the index ‘err.group_subindex(k)’. In the latter, case, the full row must be extracted for use in the histogram calculation, so the ‘:’ operator was used. The edges vector (‘x’) is shifted by half a unit to center the histogram.