TransWikia.com

(Non)LinearModelFit hangs on "trivial" linear fit (Fit succeeds)

Mathematica Asked by Effervescenza Naturale on June 29, 2021

A simplified version of my function and a simplified subset of my data (real data also has weights, if that matters, but I tested both with and without those):

myfunc[x_, a_] = 8*a*Exp[-8/Sqrt[x]]/(Sqrt[x]*Abs[WhittakerW[E*I/Sqrt[x], 
                                                             1/2, 8/10*I*Sqrt[x]*a]]^2)
mydata = {{4.5, 195}, {2.9, 175}, {2.1, 95}}

myfunc also needs a global rescaling. If I want to fit both parameter a and the scale, NonlinearModelFit works great out of the box and quickly (a Timing on a fresh kernel returned 0.6 seconds):

NonlinearModelFit[mydata, scale*myfunc[x, a], {scale, a}, {x}]

If, however, I only want to find the best scaling for an a value of my choice, this:

NonlinearModelFit[mydata, scale*myfunc[x, 4], {scale}, {x}]

takes forever in my computer (hangs until I abort). I am surprised: why the "complicated" fit works so easily, and this "trivial" linear one not?

LinearModelFit is supposedly the right tool here, but:

LinearModelFit[mydata, {myfunc[x, 4]}, {x}, IncludeConstantBasis -> False]

hangs in the same way (I let it run for about 1.5 hours, then aborted). There must be something wrong I am doing.
I tried using 4.0 instead of 4 for a, enclosing with an N[], adding WorkingPrecision -> MachinePrecision, suggesting a good starting value for scale, but I get the same result.

NonlinearModelFit works adding Method -> "NMinimize", although it is slow (Timing gave 19 seconds), and raises a warning

NonlinearModelFit::lmnl: The model (32 E^(-(8/Sqrt[x])) scale)/(Sqrt[x] Abs[WhittakerW[(I E)/Sqrt[x],1/2,(16 I Sqrt[x])/5]]^2) is linear in the parameters {scale}, but a nonlinear method or non-Euclidean norm was specified, so nonlinear methods will be used.

I do not really know what NMinimize is doing here, I guess it applies the "right" simplification to myfunc, so the computation is easier?

Finally, after wasting one day, Fit just works:

Fit[mydata, {myfunc[x, 4]}, {x}]

This reinforces my impression that I just need to apply some magical simplification, that Fit always employs, that NonlinearModelFit uses only for the two-parameter fit, and that Method -> "NMinimize" will cause. Is LinearModelFit evaluating myfunc much more times than needed? I am clueless.

LinearModelFit has several niceties that I would rather not give up for what appears to be such a stupid issue. Furthermore, I am probably going to learn a new important Mathematica quirk from this. Thus:

Why the above code, using (Non)LinearModelFit for a simple scale*f[x] model, hangs, and how should I modify it to make it work?


For completeness, I am adding some tests on my system (12.2.0.0) after @xzczd suggestions in the comments. This one:

Clear[undefinedname]
FindFit[{undefinedname}, b WhittakerW[I, 1, I], b, {x}]

hangs, regardless of whether some data is passed. LinearModelFit[{undefinedname},...] instead won’t get tricked and stop immediately. However, passing real data, those do not hang:

FindFit[mydata, b WhittakerW[I, 1, I * 1.0], ... ]
LinearModelFit[mydata, Abs[WhittakerW[I, 1, I * 1.0]], ... ]
LinearModelFit[mydata, Abs[WhittakerW[I, 1, I]], ...] (*5x slower than the previous one*)

NonlinearModelFit[mydata, scale*myfunc[x, 4], {scale}, {x},
                  Method -> "ConjugateGradient"] (*Hangs with InteriorPoint*)

It sounds like both the linear fitter and InteriorPoint get confused when they cannot explicitly evaluate WhittakerW (as it always is when the function depends on x).

One Answer

I hope that your real data consists of more than 3 data points when fitting a function with 3 parameters (a, b, and an error variance).

While I don't know why the functions hang, for your function there is a simple fix. Because myfunc[x, 4] does not change with the scaling parameter, just replace the predictor with the values of the function. Then both NonlinearModelFit and 'LinearModelFit` work fine.

myfunc[x_, a_] = 8*a*Exp[-8/Sqrt[x]]/(Sqrt[x]*
     Abs[WhittakerW[E*I/Sqrt[x], 1/2, 8/10*I*Sqrt[x]*a]]^2)

(* Replace the predictor with the result of the function *)
mydata = {{4.5, 195}, {2.9, 175}, {2.1, 95}};
mydata[[All, 1]] = myfunc[#, 4] & /@ mydata[[All, 1]];

(* Fit the model to the data *)
NonlinearModelFit[mydata, b x, {{b, 16}}, x] // Normal
(* 18.4684 x *)
LinearModelFit[mydata, x, x, IncludeConstantBasis -> False] // Normal
(* 18.4684 x *)

I know that doesn't answer the "why" part of your question.

Answered by JimB on June 29, 2021

Add your own answers!

Ask a Question

Get help from others!

© 2024 TransWikia.com. All rights reserved. Sites we Love: PCI Database, UKBizDB, Menu Kuliner, Sharing RPP