Esfit minimisation for curry - chi*T vs T and M vs H data

[thanasis]

I have been trying our curry's new capabilities for simultaneous fits to chi*T vs T and m vs H data.

While things are working apparently very well, I am not so sure as to the way the error function is calculated.
From the documentation I read:

esfit computes the error of the simulated spectrum using the root-mean-square-deviation (rmsd, i.e. the square root of the mean of the square of the deviations), where the deviation is the difference between the experimental and the simulated spectrum.

In other words, RMSD(i) = sqrt{avg[(yexp(i)-ycalc(i)]^2}, if I understand correctly.

In EPR spectra, this should be perfectly good, since absolute intensities are not important and spectra can be normalised. Even multiple spectra can be normalised and fitted simultaneously.

However, for SQUID data, I believe this could be problematic. We may observe that, from its definition, RMSD has units. And while this is not an issue for homogeneous chi*T vs T or M vs H data, when we try to fit the two types simultaneously, the two parts of the dataset will be weighed differently according to their units. E.g. chi*T in SI units can be in the range of 10e-5, whereas magnetization in CGS units can be in the 1e0 range. Clearly, the fitting function will prioritize a good fit to the M vs H data, as small relative deviations will be much larger in absolute terms. And if we switch to other combinations of units, the weighing will favor each dataset differently.

I would therefore make the suggestion that, especially for curry, the error function be calculated as [yexp(i)-ycalc(i)]^2/yexp(i)^2, or something similar, so that it becomes unitless, thereby allowing a meaningful error calculation during simultaneous fits of dissimilar data sets.

Does that make any sense?

Thanks for considering!

[katarkon]

For my experience, the fitting procedure for different datasets requires concatenation of the experimental data and writing a custom simulation function anyway. I think there are no problem to rescale both experimental and simulated data via rescale function.

[thanasis]

For EPR spectra, yes. But there is information in the absolute values of magnetic susceptibility/magnetization data that we would lose by rescaling.

[katarkon]

Yes, I'm familiar with the restriction that the magnetic data should not be automatically rescaled with esfit. I propose to rescale the initial experimental data and transfer the rescale factor to custom simulation function (which will be used anyway). So, the simulated dataset have to be rescaled by this function and the results should be correct.

[thanasis]

I agree, this is perfectly feasible and would solve this problem to a great extent.

Still, such an approach would still entail certain limitations:
-M vs H data at different temperatures can differ very widely in their saturation values (e.g. from 2 to 20 K). Thus, each dataset would have to be scaled individually.
-Within a dataset, data points with low values (e.g. magnetizations at low fields, low-T susceptibilities of AF systems or strong-zfs systems etc) would be weighed lower than data points with higher values. While this may not be such a serious issue when the data points represent spectral absorptions,* this might lead to underestimation of low-T phenomena (e.g. zfs) or low-H phenomena (e.g. steps in magnetization).

So, while computationally your approach is perfectly OK, its physical and statistical consequences would have to be checked on a per system basis. On the other hand, dividing with the experimental y value is simple, computationally undemanding, and yields a relative (and dimensionless) error which assigns an equal weight to all points of the dataset.

*Though there would be a case to be made for weak satellite features that carry information, but I let this aside for the current discussion.

[katarkon]

I think we can ask the developers to write a custom version of esfit supporting such feature. The best way is usage of additional fitting options like FitOpt.RMS='absolute'(by default) or FitOpt.RMS='relative' (RMS error is scaled according to experimental datapoint value). This feature also may be useful for more accurate fitting of EPR spectra with weak peaks (like isotope satellites or large amount of almost equivalent nuclei).
But I'm afraid, some problems may occur however.
- Termination criteria related with RMSD should be probably reconsidered
- The operability of dedicated least-squares methods like Levenberg-Marquardt is doubtful

[thanasis]

I agree!

This would definitely be a welcome addition, even if it is not applicable (at least initially) with all types of least-square fitting algorithms.

[Stefan Stoll]

As katarkon writes, if you want to fit two datasets simultaneously and have complete control over the relative scaling, a custom simulation function should work for any scenario. You have full control over the scaling of the experimental data anyway, and the custom function gives you a way to adjust the scaling of the simulated datasets. This way, the error surface can be controlled quite closely, even though the error function is hardwired into esfit. Also keep in mind to use Opt.Scaling = 'none'.

If you have a clean idea how to extend the interface of esfit to simplify usage in such cases, we are happy to consider it.

[thanasis]

With ES becoming open source, I thought of taking another look at this question.

So, I looked at esfit and found the calculation of the RMSD:

Code: Select all

Residuals = calculateResiduals(BestSpecScaled(:),FitData.ExpSpecScaled(:),FitOpts.TargetID);
Residuals = Residuals.'; % col -> row
rmsd = sqrt(mean(Residuals.^2));

What if we replaced this code with:

Code: Select all

Residuals = calculateResiduals(BestSpecScaled(:),FitData.ExpSpecScaled(:),FitOpts.TargetID);
Residuals = Residuals.'; % col -> row
if FitOpt.RMS == 'absolute'
rmsd = sqrt(mean(Residuals.^2 ./ FitData.ExpSpecScaled.^2));
else
rmsd = sqrt(mean(Residuals.^2));
end

Assuming the structure FitData.ExpSpecScaled contains the experimental data as a row, the resulting residuals should be unitless relative quantities, which would allow a rigorous fitting of any two different quantities (we could even imagine simultaneous EPR and magnetic susceptibility fits if we go crazy).
It's only a question of incorporating a FitOpt.RMS option as proposed by Katarkon.

I haven't tried to compile Easyspin, and haven't tried github (other than cloning projects), but if this is in the right direction I could try it out.

What do you think?