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
).
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
Get help from others!
Recent Answers
Recent Questions
© 2024 TransWikia.com. All rights reserved. Sites we Love: PCI Database, UKBizDB, Menu Kuliner, Sharing RPP