3
$\begingroup$

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

$\endgroup$
3
  • $\begingroup$ Looks like a bug to me. LinearModelFit doesn't hang in v9. A simpler example that hangs: Clear[aha, myfunc]; myfunc[x_?NumericQ, a_?NumericQ] := WhittakerW[I, 1, I a]; FindFit[{aha}, scale myfunc[x, 1], scale, {x}]. BTW, this sampe doesn't hang in v9.0.1, but if I change myfunc[x, 1] to myfunc[x, 3.47303], FindFit hangs. (Yeah myfunc[x, 3.47302 doesn't hang in v9. ) Have you reported this to WRI? $\endgroup$
    – xzczd
    Jan 6 at 15:48
  • $\begingroup$ This seems to be related to the method InteriorPoint. All other available methods give an output within a minute: {"Gradient", "ConjugateGradient"(*, "InteriorPoint"*), "QuasiNewton", "Newton", "NMinimize", "LevenbergMarquardt"}. $\endgroup$
    – xzczd
    Jan 6 at 16:09
  • $\begingroup$ @xzczd Very interesting. It took me a while before realizing you really meant to pass an undefined {aha}, but it does hang! I updated the question with a couple of tests. Good catch about simply changing method, although I have to give up linear fitting that way. No I haven't reported it yet, I was assuming there was some simple obscure option to pass to solve the problem. I will now. $\endgroup$ Jan 6 at 17:29
3
+50
$\begingroup$

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.

$\endgroup$
3
  • $\begingroup$ Yes, those 3 points are just for allowing people here to play around themselves if they wish, real data has more points. $\endgroup$ Jan 6 at 9:56
  • $\begingroup$ Good idea, before asking I was changing y = b f(x) to z = (y/f(x)) = b, but doing instead y = b X = b (f(x)) seems better, when adding errors to y they stay unaltered, and one can easily switch to more general linear fits. $\endgroup$ Jan 6 at 10:02
  • $\begingroup$ Hmm… I didn't expect the post is so unpopular. (My bounty doesn't even attract an upvote! ) Anyway, here's already an answer that meet the goal, so +6 :) . $\endgroup$
    – xzczd
    Jan 24 at 5:19

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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