10
$\begingroup$

I want to do a very complicated fit (it involves a product of polynomials and non-linear functions).

I read about NonlinearModelFit and thought it would be a good idea to use it. I started with a toy example:

$Version
"12.3.1 for Microsoft Windows (64-bit) (June 24, 2021)"

Clear[target, discretized, approx, min]
target[x_] := Sin[5 x] + RandomReal[{0, x}]/10
discretized = Table[{x, target[x]}, {x, 0, 1, 0.1}]
approx[x_] = Sin[a x];
min = NonlinearModelFit[discretized, approx[x], {a}, x] // Normal
Show[ListPlot[discretized], 
 Plot[min, {x, 0, 1}, PlotStyle -> {Red, Dashed}]]

I get the result -Sin[0.303389 x], which clearly is a very bad fit.

enter image description here

What happened? How can I address this? And if NonlinearModelFit can't solve these problems, what should I use instead?

Improve this question
$\endgroup$
4
  • 4
    $\begingroup$ Nonlinear optimization problems almost never just magically work without a push in the right direction. This is even especially true for problems with multiple local minima like this one. You need to give a decent starting guess for a, like: NonlinearModelFit[discretized, approx[x], {{a, 4}}, x]. Providing constraints on your fit parameters is highly recommended as well. $\endgroup$
    – Sjoerd Smit
    Commented Aug 8, 2022 at 7:04
  • $\begingroup$ @SjoerdSmit makes sense. Can you please post this as an aswer? $\endgroup$
    – infinitezero
    Commented Aug 8, 2022 at 7:13
  • 1
    $\begingroup$ @infinitezero Since the problem doesn't seem to be with the function itself, perhaps a more descriptive title might be useful for future readers? $\endgroup$
    – Hans Olo
    Commented Aug 8, 2022 at 14:10
  • $\begingroup$ @Hans Olo feel free to find a more descriptive title. $\endgroup$
    – infinitezero
    Commented Aug 8, 2022 at 15:00

3 Answers 3

Reset to default
13
$\begingroup$

Nonlinear optimization problems almost never just magically work without a push in the right direction. This is especially true for problems with multiple local minima like this one. You need to give a decent starting guess for a, like:

NonlinearModelFit[discretized, approx[x], {{a, 4}}, x].

Providing constraints on your fit parameters is highly recommended as well.

As the other answer shows, adding more freedom to your fit can also help since that freedom will prevent the fit from getting stuck in local minima. But this is not what you always want, obviously. Sometimes you want to fit a restricted function, in which case you should help the algorithm along as much as you can.

Improve this answer
$\endgroup$
1
  • $\begingroup$ Thanks for pointing this out. I somehow assumed Mathematica will do something "magical", because the example in the help section where they do a Log fit works without initial parameters and I didn't read the top section carefully enough to see that you can provide starting points. $\endgroup$
    – infinitezero
    Commented Aug 8, 2022 at 9:25
14
$\begingroup$

We need to select another fit function( shift the function Sin[a*x] to Sin[a (x + p)] + q)

approx[x_] = Sin[a (x + p)] + q;
min = NonlinearModelFit[discretized, approx[x], {a, p, q}, x] // Normal
Show[ListPlot[discretized], 
 Plot[min, {x, 0, 1}, PlotStyle -> {Red, Dashed}]]

enter image description here

Improve this answer
$\endgroup$
4
  • $\begingroup$ This baffles me. Why would fitting a more complicated model than the model that generates the data not result in a similar wrong fit as the "more appropriate" model? Clearly this works in this case. I just don't understand why. $\endgroup$
    – JimB
    Commented Aug 8, 2022 at 17:21
  • $\begingroup$ @JimB What can happen (and indeed often does) is that a simple model gets stuck in a local optimum. A more complicated model has more freedom to explore the space of possible curves, so what was previously a local optimum becomes a saddle point. $\endgroup$
    – David Zhang
    Commented Aug 10, 2022 at 7:48
  • $\begingroup$ There is a real art to introducing artificial "slack" parameters into nonlinear optimization problems to provide just the right amount of extra freedom without turning the problem into an under-constrained mess. Unfortunately, I don't know of any good, principled ways to know how much slack to introduce ahead of time; you have to experiment with your particular problem. (Also note that this is an unrelated meaning of "slack" to "slack variables" in linear programming.) $\endgroup$
    – David Zhang
    Commented Aug 10, 2022 at 7:51
  • 1
    $\begingroup$ @DavidZhang That makes sense. But then wouldn't one need to go back to the simpler model and use the estimate of $a$ as a starting value? Also, in this case both other the "extraneous/slack" parameters are statistically significant. That also is not comforting (although the estimators for $a$ and $p$ are nearly perfectly correlated). In short, I see that adding the slack parameters gets the algorithm "unstuck". But one just can't stop there. That would mean the "numerical instability" is highly (and inappropriately) influencing inferences about the parameters. $\endgroup$
    – JimB
    Commented Aug 10, 2022 at 16:19
9
$\begingroup$

There are a few issues with the fitting. Some issues are issues with fitting in general (i.e., no matter what software package is used) and some issues are caused by you. First the issues caused by you:

  1. Maybe it's different in your field but I find it rare that errors have a uniform distribution.
  2. You have the errors increase with the predictor value (which is fine) but you then attempt a fit with NonlinearModelFit that assumes the variance of the errors is constant for all values of the predictors.

I understand that this is a toy example but understand you have put in issues that can cause a lack of fit that are not the fault of NonlinearModelFit.

The function being fit (a sine curve) typically requires good starting values (as NonlinearModelFit uses an iterative fitting procedure that attempts to find parameters that minimize a sum of squares). The sum of squares being minimized for this example fluctuates wildly (i.e. not in a nice convex manner):

(* Generate data *)
SeedRandom[12345];
target[x_] := Sin[5 x] + RandomReal[{0, x}]/10
discretized = Table[{x, target[x]}, {x, 0, 1, 0.1}]
(* Calculate sum of squares for various values of the parameter a *)
sumOfSquares = (discretized[[All, 2]] - 
     Sin[a discretized[[All, 1]]])^2 // Total
(* 0. + (0.482725 - Sin[0.1 a])^2 + (0.857126 - Sin[0.2 a])^2 + (1.0104 -
    Sin[0.3 a])^2 + (0.918241 - Sin[0.4 a])^2 + (0.621625 - 
   Sin[0.5 a])^2 + (0.185401 - Sin[0.6 a])^2 + (-0.30125 - 
   Sin[0.7 a])^2 + (-0.69353 - Sin[0.8 a])^2 + (-0.968016 - 
   Sin[0.9 a])^2 + (-0.92088 - Sin[1. a])^2 *)

(* Plot sum of squares vs value of a *)
Plot[sumOfSquares, {a, -20, 20}, PlotRange -> {All, {0, Automatic}}, Frame -> True, 
 FrameLabel -> (Style[#, Bold, 18] &) /@ {"a", "Sum of squares"}]

Sum of squares vs value of a

From the above plot it appears that getting an estimate of $a$ around 5 minimizes the sum of squares (which clearly matches the value used to generate the data). However, because there are many local minima (as pointed out by @SjoerdSmit), the algorithm can stop at the wrong value.

If we look at starting values from -23 to 20, there are 7 (essentially) unique results. Here is a plot of the starting value vs. the resulting estimate of $a$:

results = 
  Table[{a0, a, sumOfSquares[a], sumOfSquares[a0]} /. 
    NonlinearModelFit[discretized, Sin[a x], {{a, a0}}, x][BestFitParameters"],
  {a0, -23, 20, 1/10}];
(* Round the estimates of a to account for not stopping at the exact same value *)
results[[All, 2]] = Round[results[[All, 2]], 0.001];

(* Separate results by the unique values of the estimates of a *)
uniqueA = DeleteDuplicates[results[[All, 2]]];
r = Table[Select[results, #[[2]] == uniqueA[[i]] &], {i, Length[uniqueA]}];

(* Plot results *)
Show[ListPlot[r[[All, All, {1, 4}]], Joined -> True, 
  PlotRange -> {All, {0, Automatic}}, Frame -> True, 
  FrameLabel -> (Style[#, Bold, 18] &) /@ {"a", "Sum of squares"}],
 ListPlot[r[[All, All, {2, 3}]], PlotStyle -> PointSize[0.02]], 
 PlotRangeClipping -> False]

Results that depend on the starting value

So for this particular model and data generation process, starting values for $a$ should be between 2 and 8.

I'll also note that your not-completely described real model ("a product of polynomials and non-linear functions") might also be a linear model if that model is "linear in the parameters" even though there are nonlinear functions of the predictors. In that case LinearModelFit would be more stable (as it doesn't require - or even allow - starting values).

Improve this answer
$\endgroup$
3
  • $\begingroup$ Thanks for your insights and of course you're completely right! As mentioned in another comment, I got carried away by the mathematica example that does not require starting values and my brain just went blank and stopped thinking about it. As a side note: my actual fit parameters appear in terms $\propto 1/x^2$ so it definitely is non-linear. $\endgroup$
    – infinitezero
    Commented Aug 8, 2022 at 18:02
  • $\begingroup$ $a_0+a_1 \log(x)+a_2/x^2$ is a linear function that can be fit with LinearModelFit. Or do you mean something different by "parameters appear in terms of $\propto 1/x^2$"? $\endgroup$
    – JimB
    Commented Aug 8, 2022 at 18:06
  • $\begingroup$ Sorry for being ambigious. I meant something like $a_0 \frac{c}{d^2+(t-b_0)^2} f_0(x)$ where $a_i$ and $b_i$ are element of my parameter set and $d$ and $t$ are constants. $\endgroup$
    – infinitezero
    Commented Aug 8, 2022 at 18:09

Your Answer

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

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