2
$\begingroup$

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}}

And the "Bad" data:

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?

Improve this question
$\endgroup$
2
  • $\begingroup$ @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. $\endgroup$
    – JimB
    Jan 3, 2019 at 20:07
  • $\begingroup$ Show a working example of how this looks. $\endgroup$
    – Alex Trounev
    Jan 3, 2019 at 21:05

5 Answers 5

Reset to default
5
$\begingroup$

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" ] 
Show[{ListPlot[BadFitDataU], 
Plot[Normal[bad], {\[Nu], Min[BadFitDataU [[All, 1]]],Max[BadFitDataU [[All, 1]]]} , PlotRange -> All ]}]

enter image description here

detects a peak!

Improve this answer
$\endgroup$
5
$\begingroup$

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

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

ClearAll["Global`*"]

{m1, m2} = Mean[BadFitData];
{\[Sigma]1, \[Sigma]2} = StandardDeviation[BadFitData];
data = BadFitDataU = 
   Transpose[{(BadFitData[[All, 1]] - 
        m1)/\[Sigma]1, (BadFitData[[All, 2]] - m2)/\[Sigma]2}];

  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]

enter image description here

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..

Improve this answer
$\endgroup$
4
$\begingroup$

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.

enter image description here

Improve this answer
$\endgroup$
6
  • 2
    $\begingroup$ 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. $\endgroup$
    – JimB
    Jan 4, 2019 at 1:18
  • 2
    $\begingroup$ 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). $\endgroup$
    – JimB
    Jan 4, 2019 at 1:25
  • $\begingroup$ 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... :) $\endgroup$
    – MikeY
    Jan 4, 2019 at 4:43
  • 1
    $\begingroup$ Not if you're doing medical research or building bridges or weapons. You're entitled to your artistry but not to bogus statistical theory. $\endgroup$
    – JimB
    Jan 4, 2019 at 5:32
  • $\begingroup$ It’s OK, I’m a Bayesian and I had a prior. $\endgroup$
    – MikeY
    Jan 4, 2019 at 12:25
3
$\begingroup$

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]}}]]

Good data fit

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]}}]]

Bad data fit

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 *)
Improve this answer
$\endgroup$
2
  • $\begingroup$ 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. $\endgroup$
    – MikeY
    Jan 4, 2019 at 19:23
  • 1
    $\begingroup$ @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). $\endgroup$
    – JimB
    Jan 4, 2019 at 19:46
2
$\begingroup$

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/\
MathematicaForPrediction/master/MonadicProgramming/\
MonadicQuantileRegression.m"]

qrObj =
  QRMonUnit[BadFitData]⟹
   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];

enter image description here

Improve this answer
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

Not the answer you're looking for? Browse other questions tagged or ask your own question.