# Making batch data fitting robust — why does NonlinearModelFit fail occasionally?

I have over 1000 spectra I would like to batch fit and then use the results for further study.

I'm not fitting anything especially complicated just a Lorentzian of the form, $$y(x) = A \frac{\gamma^{2} /2}{(x - x_{0})^{2} + \gamma^{2}/4} + y_{0}.$$

The problem I have is that some of my spectra fit quite happily and with good/reasonable results. However some spectra despite looking very similar and the respective fit being initialised in a reasonable way produce very odd fit results. Here is the "Good" data:

GoodFitData = {{29.642582062500000000000,
4.3190002032711183*10^-7}, {29.6425820937500000000000,
5.3406344114652258*10^-7}, {29.642582125000000000000,
4.3050696173846918*10^-7}, {29.642582156250000000000,
3.8476329392273719*10^-7}, {29.642582187500000000000,
3.0946371958194098*10^-7}, {29.642582218750000000000,
3.3232729264197710*10^-7}, {29.642582250000000000000,
6.3905435090395855*10^-7}, {29.6425822812500000000000,
8.753370135668026*10^-7}, {29.642582312500000000000,
8.3216510981331772*10^-7}, {29.642582343750000000000,
5.6798059206608427*10^-7}, {29.642582375000000000000,
3.2057831588605133*10^-7}, {29.642582406250000000000,
6.0006018141506731*10^-7}, {29.642582437500000000000,
5.9220334443374970*10^-7}, {29.642582468750000000000,
5.9043428760173520*10^-7}, {29.6425825000000000000000,
5.3142937509568503*10^-7}, {29.6425825312500000000000,
1.1512362985808565*10^-6}, {29.642582562500000000000,
2.2005730090913985*10^-6}, {29.642582593750000000000,
1.3061501910622896*10^-6}, {29.642582625000000000000,
7.0306201913894104*10^-7}, {29.642582656250000000000,
7.5940811728371987*10^-7}, {29.642582687500000000000,
8.896327027735458*10^-7}, {29.6425827187500000000000,
5.1906915196157101*10^-7}, {29.642582750000000000000,
5.0793870783619110*10^-7}, {29.642582781250000000000,
3.9490581635532259*10^-7}, {29.642582812500000000000,
3.2285197574563654*10^-7}, {29.642582843750000000000,
5.0979050094687087*10^-7}, {29.642582875000000000000,
5.7917798046612257*10^-7}, {29.642582906250000000000,
2.9746088835309100*10^-7}, {29.642582937500000000000,
2.8446466006038211*10^-7}, {29.642582968750000000000,
5.1303286544108108*10^-7}, {29.642583000000000000000,
4.4818646841234343*10^-7}, {29.6425830312500000000000,
3.3482541855830208*10^-7}, {29.6425830625000000000000,
2.2472479916643503*10^-7}}


BadFitData = {{29.642582375000000000000,
1.4966773070074960*10^-7}, {29.6425823906250000000000,
3.2796187771352530*10^-7}, {29.642582406250000000000,
3.3407923869075059*10^-7}, {29.6425824218750000000000,
1.4881722372170584*10^-7}, {29.642582437500000000000,
1.5636813428753172*10^-7}, {29.642582453125000000000,
4.1457809704998334*10^-7}, {29.642582468750000000000,
5.1785323732788303*10^-7}, {29.642582484375000000000,
4.0880319580912151*10^-7}, {29.6425825000000000000000,
2.4201807866572215*10^-7}, {29.6425825156250000000000,
2.8251493555106467*10^-7}, {29.6425825312500000000000,
4.8387544256366870*10^-7}, {29.642582546875000000000,
5.5329474183042744*10^-7}, {29.642582562500000000000,
6.5744418237327978*10^-7}, {29.642582578125000000000,
7.9552641010675053*10^-7}, {29.642582593750000000000,
4.4648275313948193*10^-7}, {29.642582609375000000000,
5.9516896108447851*10^-7}, {29.642582625000000000000,
2.4332278859121748*10^-6}, {29.642582640625000000000,
1.9634656151700146*10^-6}, {29.642582656250000000000,
1.7604070459744035*10^-7}, {29.642582671875000000000,
2.2592814520073882*10^-7}, {29.642582687500000000000,
2.6774361092830216*10^-7}, {29.6425827031250000000000,
3.3119779204170161*10^-7}, {29.6425827187500000000000,
5.1003549704521799*10^-7}, {29.6425827343750000000000,
2.7701066539193851*10^-7}, {29.642582750000000000000,
2.5924933518664618*10^-7}, {29.642582765625000000000,
4.6290702798493546*10^-7}, {29.642582781250000000000,
4.5979200634303145*10^-7}, {29.642582796875000000000,
6.4157516355120082*10^-7}, {29.642582812500000000000,
4.2154540563855831*10^-7}, {29.642582828125000000000,
2.3994878809517319*10^-7}, {29.642582843750000000000,
4.8672163908576681*10^-7}, {29.642582859375000000000,
5.0605896869103541*10^-7}, {29.642582875000000000000,
4.4161357956856321*10^-7}}


My fit code is:

LorentzFunction[A_,\[Gamma]_,y0_,\[Nu]0_,\[Nu]_] =  A (\[Gamma]^2/2)/((\[Nu] - \[Nu]0)^2 + \[Gamma]^2/4) + y0;
NonlinearModelFit[
CutPeakData,
LorentzFunction[A,\[Gamma],y0,\[Nu]0,\[Nu]],
{{A,PeakSignal},{\[Gamma],2*10^-8},{\[Nu]0,PeakFrequency},{y0,5*10^-7}},\[Nu],
MaxIterations->1000
];


I initialise the fit parameters by extracting the peak $$y$$ and $$x$$ values for A and \[Nu]0 respectively. For \[Gamma] I set 2*10^-8 and y0 I take the mean of the signal floor both sides of the peak.

I've tried various approaches such as the ones described here -- which is a really excellent guide and has worked for me with other data sets (exponential fits) -- and playing with the number of MaxIterations. None seem to help the BadFitData despite the data itself not looking that bad, at least by eye -- thoughts?

• @Feyre 's answer (mathematica.stackexchange.com/questions/139038/…) is excellent. (And I added a couple suggestions, too). Also, to note that this isn't just a NonlinearModelFit or Mathematica issue. SAS, R, SPSS, etc. all have similar issues. It's just the nature of iterative fitting AND users not following guidelines as suggested in @Feyre's answer. – JimB Jan 3 at 20:07
• Show a working example of how this looks. – Alex Trounev Jan 3 at 21:05

It seems to be a "scaling" problem(I only consider the bad data):

BadFitDataU = Standardize[BadFitData] (*standardization: mean=0, standarddeviation=1*)


Fitting this data using Method->"NMinimize" without any additional assumptions

bad = NonlinearModelFit[BadFitDataU,LorentzFunction[A, \[Gamma],y0, \[Nu]0, \[Nu]], {A, \[Gamma], \[Nu]0 , y0 }, \[Nu],Method -> "NMinimize" ]


detects a peak!

Here is simpler model. Let's use @Ulrich;s scaled data..

$$f(x)=\frac{a}{(x-\text{x0})^2+b}+c$$

ClearAll["Global*"]

f[x_] := a/((x - x0)^2 + b) + c;

n = Length@data;
parameters = {a, b, c, x0};
cost = \!$$\*UnderoverscriptBox[\(\[Sum]$$, $$i = 1$$, $$n$$]
\*SuperscriptBox[$$(f@data[[i, 1]] - data[[i, 2]])$$, $$2$$]\) ;(* Cost function to Minimize *)
fit = NMinimize[cost, parameters, Method -> "DifferentialEvolution",
MaxIterations -> 5000];

Thread[{a, b, c, x0} = parameters /. Last@fit];

Plot[f[x], {x, data[[1, 1]], data[[-1, 1]]}, Epilog :> Point[data],
PlotStyle -> {Orange, Thick}, PlotRange -> All, Frame -> True,
Axes -> False]


Of course you can use NonlinearModelFit

nlm = NonlinearModelFit[data, f[x], {a, b, c, x0}, x,
Method -> "NMinimize"]

Plot[nlm // Normal, {x, data[[1, 1]], data[[-1, 1]]},
Epilog :> Point[data], PlotStyle -> {Orange, Thick},
PlotRange -> All, Frame -> True, Axes -> False]


Same picture..

Scaling is critical in my experience.

Something else I've done when in extremis to kick-start the nonlinear model fitting/debugging process is just augment or smooth the data. NonLinearModelFit likes lots of points. Try...

CutPeakData = Join[ BadFitData , MovingAverage[BadFitData,2] ]


...which just adds points between the points. Then run the fit. The graphic shows the original data points, the fit is to the augmented data. It can allow you to at least see if your fit is shaped right. Super-easy to do too. Can use the fitted model parameters as new start points.

• Adding in smoothed data as if it were real data is at best inappropriate. I know it happens but there is no theoretical reason for it. In fact, what you get at best is overly optimistic estimates of precision. (Although the folks that I know that do this regularly don't seem to care about estimates of precision in any event.) But I agree that scaling is critical. – JimB Jan 4 at 1:18
• I've retracted my -1 vote. I see you've just considered using that blasphemous technique just to get starting values. (I'm not convinced that's always helpful but at least you've mentioned that one needs to do the final fit with the real data). – JimB Jan 4 at 1:25
• JimB, plenty of room for art in the science of fitting models. Sometimes you have to just let theory hold your beer while you get to engineering. Whatever gets you off top dead center... :) – MikeY Jan 4 at 4:43
• Not if you're doing medical research or building bridges or weapons. You're entitled to your artistry but not to bogus statistical theory. – JimB Jan 4 at 5:32
• It’s OK, I’m a Bayesian and I had a prior. – MikeY Jan 4 at 12:25

I think we're all in agreement that standardizing/scaling of either the response variable or the predictor variable or both is warranted. (Sometimes scaling the parameters is also necessary.)

What's different about this answer is that standardization is by subtracting the minimum and dividing by the range to get all values between 0 and 1. Also, I show that doing so doesn't change the function being fit and how to get back to the original parameters. (The point here is that you just can't always just replace the original values with the standardized values into the function and fit the same model.) Also, I construct initial estimates for all parameters from the data.

Two of many ways of standardizing are by (1) subtracting the mean and dividing by the standard deviation (roughly getting most of the values between -4 and 4), and (2) subtracting the minimum and dividing by the range (getting all of the values between 0 and 1).

I'll use (2) here but (1) will work fine, too.

fit[data_] := Module[{νmin, νmax, νrange, ymin, ymax, yrange, stdData},

(* Standardize response and predictor *)
{νmin, νmax} = MinMax[data[[All, 1]]];
{ymin, ymax} = MinMax[data[[All, 2]]];
νrange = νmax - νmin;
yrange = ymax - ymin;
stdData = Transpose[{(data[[All, 1]] - νmin)/νrange, (data[[All, 2]] - ymin)/yrange}];

(* Get starting values based on the data (the "s" attached to the parameters and initial values
is to indicate "standardized" *)
y0sInit = Min[stdData[[All, 2]]];
asInit = (Max[stdData[[All, 2]]] - y0sInit)/2;
ν0sInit = Select[stdData, #[[2]] == Max[stdData[[All, 2]]] &][[1, 1]];
midv = (Min[stdData[[All, 1]]] + ν0sInit)/2;
midy = Select[Sort[stdData, #1[[1]] < #2[[1]] &], #[[1]] >= midv &][[1, 2]];
γsInit = Abs[γ /. Solve[LorentzFunction[asInit, γ, y0sInit, ν0sInit,  midv] == midy, γ][[1]]];

(* Find maximum likelihood estimates *)
nlm = NonlinearModelFit[stdData, LorentzFunction[as, γs, y0s, ν0s, ν],
{{as, asInit}, {γs, γsInit}, {ν0s, ν0sInit}, {y0s, y0sInit}}, ν, MaxIterations -> 50000];

(* Return results *)
{nlm, {A -> as yrange, γ -> γs νrange, y0 -> ymin + y0s yrange, ν0 -> νmin + ν0s νrange} /.
nlm["BestFitParameters"], {νmin, νmax, νrange, ymin, ymax, yrange}}]


Now to try it on the two datasets:

f = fit[GoodFitData];
{νmin, νmax, νrange, ymin, ymax, yrange} = f[[3]];
Show[ListPlot[GoodFitData, PlotRange -> All, ImageSize -> Large],
Plot[ymin + yrange f[[1]][(ν - νmin)/νrange], {ν, νmin, νmax}, PlotPoints -> 100,
PlotRange -> {All, {ymin, ymin + yrange f[[1]][((ν0 /. f[[2]]) - νmin)/νrange]}}]]


f = fit[BadFitData];
{νmin, νmax, νrange, ymin, ymax, yrange} = f[[3]];
Show[ListPlot[BadFitData, PlotRange -> All, ImageSize -> Large],
Plot[ymin + yrange f[[1]][(ν - νmin)/νrange], {ν, νmin, νmax}, PlotPoints -> 100,
PlotRange -> {All, {ymin, ymin + yrange f[[1]][((ν0 /. f[[2]]) - νmin)/νrange]}}]]


I would claim that both datasets are problematic even without the scaling issues. The "peak" of the curves is essentially fit with just two or three points near the peak (assuming the peak actually exists). Yes, I know I'm more pessimistic about examining fits than several folks around here but with lots of noise observed in the data, relying on just one or two data points is living dangerously. Finding confidence intervals for the locations of the peaks is warranted - but that's another question.

Now...does standardization change the function being fit? The question is important because sometimes just plugging in the standardized values changes the function being fit. Here's one way to be convincing that the same function is being fit for this particular model and data:

We should have the backtransformed fit with the standardized data be the same as if we had the original fit:

LorentzFunction[A_, γ_, y0_, ν0_, ν_] =  A (γ^2/2)/((ν - ν0)^2 + γ^2/4) + y0;
standardizedModelBacktransformed = FullSimplify[ymin +
yrange LorentzFunction[A/yrange, γ/νrange, (y0 - ymin)/yrange, (ν0 - νmin)/νrange, (ν - νmin)/νrange]]
FullSimplify[standardizedModelBacktransformed - LorentzFunction[A, γ, y0, ν0, ν]]
(* 0 *)

• MMA has a Standardize[ ] function that may handle all of your standardizing out of the box. I tried it on this problem, but it got fiddly coming up with a γ start parameter. – MikeY Jan 4 at 19:23
• @MikeY Thanks. I needed to calculate the min and max for backtransforming, but that function is handy (especially because it will do all columns simultaneously). – JimB Jan 4 at 19:46

As others said in previous answers we do/experiment with data rescaling. I want to add that the choice of function basis is also important.

Below is an example using (i) Quantile Regression with B-splines and Chebyshev polynomials and (ii) Linear Regression with Chebyshev polynomials.

Import["https://raw.githubusercontent.com/antononcube/\

qrObj =
QRMonEcho[Style["Original 'bad' data summary:", Bold, Purple, Larger]]⟹
QRMonEchoDataSummary⟹
QRMonEcho[Style["Original 'bad' data,\nQuantile Regression with B-splines:", Bold, Purple, Larger]]⟹
QRMonQuantileRegression[12]⟹
QRMonPlot[PlotRange -> All]⟹
QRMonEcho[Style["Original 'bad' data,\nQuantile Regression and Linear regression with Chebyshev polynomials:", Bold, Purple, Larger]]⟹
QRMonQuantileRegressionFit[12]⟹
QRMonFit[8]⟹
QRMonPlot[PlotRange -> All]⟹
QRMonEcho[Style["Rescaled 'bad' data summary:", Bold, Purple, Larger]]⟹
QRMonRescale[]⟹
QRMonEchoDataSummary⟹
QRMonEcho[Style["Rescaled 'bad' data,\nQuantile Regression with B-splines,\nLinear regression with Chebyshev polynomials:", Bold, Purple, Larger]]⟹
QRMonQuantileRegressionFit[12]⟹
QRMonFit[8]⟹
QRMonPlot[PlotRange -> All];
`