# Sensitivity of NonlinearModelFit to model

I have come across a circumstance where NonlinearModelFit is very sensitive to the model used. I am aware that NonlinearModelFit is very dependent on the initial estimates and this dictated my choice of model -- I thought I had chosen a good model. I would like to hear comments on why my choice is poor.

I am fitting data that is a cosine wave. The two choices of model I considered are

m1 = a Cos[2 π f t] + b Sin[2 π f t];
m2 = a Cos[2 π f t + ϕ];


The first model looks better because it has one nonlinear parameter, the frequency f, while the second has frequency and phase angle ϕ. I was hoping that I could just guess the frequency and not supply estimates for a and b because they are linear.

To test these two models I used the following data based on measured values.

data = With[{a = 43.45582489316203,
f = 94.92003941300389,
ϕ = 431.155471523826},
SeedRandom[1234];
Table[{t,
a Cos[2 π f t + ϕ] + RandomReal[{-0.1, 0.1}]},
{t,13.439999656460714, 13.479799655455281, 0.0002}]
];


Here is the first fit

fit1 = NonlinearModelFit[data, m1, {{f, 100}, {a, -22}, {b, 35}}, t];
fit1["ParameterConfidenceIntervalTable"]
Show[ListPlot[data], Plot[fit1[t], {t, data[[1, 1]], data[[-1, 1]]}]]


The error

Failed to converge to the requested accuracy or precision within 100 iterations.

is produced. The standard errors are poor:

Now consider the second model

fit2 = NonlinearModelFit[data, m2, {{f, 100}, {a, 40}, {ϕ, 0.7}},t];
fit2["ParameterConfidenceIntervalTable"]


This model does converge and is a good fit although the phase is several multiplies of Pi. On a minor point changing the phase to say 3.9 results in almost the same values. Is the numerical evaluation of the trig functions an issue?

fit2 = NonlinearModelFit[data, m2, {{f, 100}, {a, 40}, {ϕ, 3.9}},
t];
fit2["ParameterConfidenceIntervalTable"]
Show[ListPlot[data], Plot[fit2[t], {t, data[[1, 1]], data[[-1, 1]]}]]


I wondered if my assumption was wrong and if there was more than one minimum for the first model. I therefore generated the error on the assumption that given an estimate of frequency the problem is a linear one and a and b can be solved using LeastSquares. This module generates the mean square error given a value of frequency.

ClearAll[err];
err[data_, f_] := Module[{tt, d, mat, a, b, fit},
tt = data[[All, 1]];
d = data[[All, 2]];
mat = {Cos[2 π f #], Sin[2 π f #]} & /@ tt;
{a, b} = LeastSquares[mat, d];
fit = a Cos[2 π f #] + b Sin[2 π f #] & /@ tt;
{f, (d - fit).(d - fit), {a, b}}
]

e1 = Table[err[data, f], {f, 40, 150, 1}];
ListPlot[e1[[All, {1, 2}]]]


As expected there is a good minimum around the correct frequency with a reasonable guessing range for just the frequency. This reinforces my idea that model 1 should be better.

What's wrong with my intuition? Why is model 2 better than model 1?

• There is nothing wrong with your first model: it just needs to do more iterations to get to convergence, as it says in the error message. Try adding MaxIterations -> 1000 as an option to NonlinearModelFit. Doing so generated the following best fit parameters for me: {f -> 94.9197, a -> -30.7143, b -> 30.7426}, with errors on the parameters: {0.00207457, 5.39375, 5.38866}. This seems in good agreement at least with the frequency you found. – MarcoB Aug 28 '15 at 16:07
• @MarcoB Thanks but I000 iterations did not do much for me- I had tried this before and still had the error. There was some improvement with 10000 iterations but the standard error was poor compared to the second model. I a using version 10.0.2 with Windows(64-bit). – Hugh Aug 28 '15 at 16:59
• I am using MMA 10.2 on Win7-64. The fitting routines must have been very recently improved then. – MarcoB Aug 28 '15 at 17:03
• For both models "LevenbergMarquardt" method is used. So the method was not the cause of difference of the model success. "LevenbergMarquardt" does only use "TrustRegion" as step control. Changing the setting on the step control did not improve convergence of the first model. The conclusion for me is that it is intuition that misleads, as both models are handeled with the same solve method. – Eisbär Sep 1 '15 at 20:22
• One of the issues likely causing the instability and the need to increase the number of iterations (as stated by @MarcoB) is that the model being fit is maybe too complex for the data. This is evident from examining the parameter correlation matrix: fit1["CorrelationMatrix"] // MatrixForm results in $\left( \begin{array}{ccc} 1. & -0.999999 & -0.999999 \\ -0.999999 & 1. & 0.999999 \\ -0.999999 & 0.999999 & 1. \\ \end{array} \right)$. – JimB Aug 23 '20 at 1:25

Intuition is sometimes tricky on fitting procedures. This is of course not a Mathematica issue, but a problem of fitting in general. You can see the problem in parameter space (hence it depends on the details of parameter space).

Defining for the residuals (square root)

Res1[ff_, aa_, bb_] := Norm[data[[All, 2]] - (m1 /. {f -> ff, a -> aa, b -> bb, t -> #} & /@data[[All, 1]])]


and plotting

GraphicsGrid[{{
Plot3D[Res1[100.1, aa, bb], {aa,-50,50}, {bb,-50,50}, MeshFunctions -> {#3 &}],
Plot3D[Res1[100., aa, bb], {aa,-50,50}, {bb,-50,50}, MeshFunctions -> {#3 &}],
Plot3D[Res1[99.9, aa, bb], {aa,-50,50}, {bb,-50,50}, MeshFunctions -> {#3 &}]
}}]


you see that the gradient in the $$(a,b)$$ projection of the parameter space complete changes direction upon small changes in frequency. On the other hand with

Res2[ff_, aa_, ϕϕ_] := Norm[data[[All,2]] - (m2 /. {f -> ff, a -> aa, ϕ -> ϕϕ, t -> #} & /@ data[[All, 1]])]


and plotting

GraphicsGrid[{{
Plot3D[Res2[100.1, aa, fi], {aa,-50,50}, {fi,-Pi,Pi}, MeshFunctions -> {#3 &}],
Plot3D[Res2[100.0, aa, fi], {aa,-50,50}, {fi,-Pi,Pi}, MeshFunctions -> {#3 &}],
Plot3D[Res2[99.9, aa, fi], {aa,-50,50}, {fi,-Pi,Pi}, MeshFunctions -> {#3 &}]
}}]


is more 1 dimensional. So you are not running in circles. While not a complete answer, I hope this gives an idea.

A note at the end. My general advice is: whenever possible redefine your model such that all parameters are on the same order of magnitude.

First Update

Concerning op's concern: The plots for model 1 look nice and quadratic (as I suggested in the second part of my question). The plots for model 2 are wild and could easily take you off in the wrong direction.

I agree, but this is only in a 2D cut of the 3D problem. Moreover, phi is restricted to mod $$2 \pi$$ Sure, there are saddle points and they actually take you off, resulting in the large phase in the end, while $$431 \mod 2\pi$$ makes $$3.9$$ a good guess. Furthermore, if you jump in the next minimum of the phase and make a phase shift of $$\pi$$, the cut in amplitude is parabolic, giving you very fast the amplitude with opposite sign.

In detail you can see what I mean If you look how Mathematica travels through your parameter space (at the moment I only have Version 6 at hand)

{fit3, steps3} =
Reap[FindFit[data, m1, {{f, 100}, {a, 8}, {b, 41}}, t,
MaxIterations -> 1000, StepMonitor :> Sow[{f, a, b}]]];
Show[Graphics3D[
Table[{Hue[.66 (i - 1)/(Length[First@steps3] - 1)],
AbsolutePointSize[7], Point[(First@steps3)[[i]]],
Line[Take[First@steps3, {i, i + 1}]]}, {i, 1,
Length[First@steps3] - 1}], Boxed -> True, Axes -> True],
BoxRatios -> {1, 1, 1}, AxesLabel -> {"f", "a", "b"}]


Here you see what I mean with going in circles. Even after $$1000$$ Iterations you are not even close as the $$(a,b)$$-minimum changes position with changes in frequency in such an unfortunate way.

If you look on the other hand at the second model you get:

{fit2, steps2} =
Reap[FindFit[data, m2, {{f, 100}, {a, 40}, {ϕ, 3.9}}, t,
StepMonitor :> Sow[{f, a, ϕ}]]];
Show[Graphics3D[
Table[{Hue[.66 (i - 1)/(Length[First@steps2] - 1)],
AbsolutePointSize[7], Point[(First@steps2)[[i]]],
Line[Take[First@steps2, {i, i + 1}]]}, {i, 1,
Length[First@steps2] - 1}], Boxed -> True, Axes -> True],
BoxRatios -> {1, 1, 1}, AxesLabel -> {"f", "a", "ϕ"}]


where it finds the amplitude quite fast, reducing the problem to 2D in phase and frequency.

Second Update

Concernings the op's question if the final result is a quadratic well.

Let us just plot the three cuts in parameter space.

{faPlot = ContourPlot[Res2[freq, amp, 434.3256],
{freq, 94.9197 - .01, 94.9197 + .01}, {amp,-43.4566-10, -43.4566 +10}],
fpPlot = ContourPlot[Res2[freq, -43.4566, phase],
{freq, 94.9197-.01,94.9197 + .01}, {phase, 434.3256-1.5,434.3256+1.5}],
apPlot = ContourPlot[Res2[94.9197, amp,phase],
{amp, -43.4566-15, -43.4566+15}, {phase,434.3256-1.5,434.3256+1.5}]}


This looks promising except for the middle graph. After a coordinate transformation, however, we get

β = 84.57;
ContourPlot[Res2[94.9197 + fff + ppp/β, -43.4566, 434.3256 - β fff + ppp],
{fff, -8.5, +8.5}, {ppp, -1.5, +1.5}]


which gives

So this looks OK as well. All is good.

Making the troublesome fit work

On StackOverflow I came across answers from Jean Jaquelin providing methods to turn non-linear fits in actual linear fits. Some information can be found here.

The point is that when looking at $$y = a \cos( \omega t) + b \sin(\omega t)$$ we know that the second derivative is $$y'' = - \omega^2 y$$. Numerical derivatives are very often critical though. Slightly better is to look at the double integration $$\int\int y = -y/\omega^2 + c t + d$$. The Integration can be performed rather easy with

cumint[ indata_ ] := Module[
{
p = Interpolation[indata] ,
timedata,
signaldata,
int
},
timedata = indata[[All, 1]];
signaldata = indata[[All, 2]];
int = Join[{0}, Table[
NIntegrate[ p[t], {t, timedata[[i]], timedata[[i + 1]]} ],
{i, 1, Length[ timedata ] - 1 } ]
];
Return[ Transpose[{ timedata, Accumulate[int] } ] ]
]


(This is my quick and dirty solution while in python using cumtrapz) This leaves us with a linear optimization for $$1/\omega^2$$, $$c$$ and $$d$$, while we are only interested in the first one. We then have

dT = Transpose[data];
tList = dT[[1]];
sList = dT[[2]];
y1 = cumint[data];
y2 = cumint[y1];
SSList = y2[[All, 2]];
GraphicsArray[{{ListPlot[ data, Joined -> True],
ListPlot[ y1, Joined -> True], ListPlot[ y2, Joined -> True]}}]


VT = {sList, tList, Table[1, Length[data]]};
V = Transpose[ VT ];
A = VT.V;
SV = VT.SSList;
AI = Inverse[ A ];
\[Alpha] = AI.SV;
w0 = Sqrt[-1/\[Alpha][[1]]];
f0 = w0/2/Pi


which gives f0 = 94.9134. With this knowledge one can make a linear fit on a and b, namely

sv = Sin[ w0 tList];
cv = Cos[w0 tList];
WT = {cv, sv};
W = Transpose[WT];
B = WT.W;
BI = Inverse[B];
SY = WT.sList;
sol = BI.SY


providing {-10.9575, 42.0522}, and

ListPlot[
{
sList, (sol[[1]] cv +  sol[[2]] sv)
}, Frame -> True, Joined -> {False, True}
]


This looks already very good. Now let's try to use this results as start parameters for the non-linear fit.

{fit3, steps3} =
Reap[FindFit[data, m1, {{f, f0}, {a, sol[[1]]}, {b, sol[[2]]}}, t,
MaxIterations -> 1000, StepMonitor :> Sow[{f, a, b}]]];
Show[Graphics3D[
Table[{Hue[.66 (i - 1)/(Length[First@steps3] - 1)],
AbsolutePointSize[7], Point[(First@steps3)[[i]]],
Line[Take[First@steps3, {i, i + 1}]]}, {i, 1,
Length[First@steps3] - 1}], Boxed -> True, Axes -> True],
BoxRatios -> {1, 1, 1}, AxesLabel -> {"f", "a", "b"}]
fit3


{f -> 94.9197, a -> -30.7143, b -> 30.7426}

Now it works. The slight modification in the frequency, however, resulted again in a quite dramatic change of amplitudes. Does it really fit? How does it look?

ListPlot[
{
data,
Table[{t, a Cos[ 2 Pi f t] + b  Sin[2 Pi f t]}, {t, data[[1, 1]],
data[[-1, 1]], 0.0001}] /. fit3
}, Joined -> {False, True}
]


It does fit and looks good. In this simple case the pure linear approach probably would have been enough. In case of noisy data it still might work, but definitively gives a good set of starting values. One also can use this results to calculate better starting values for the solution using phases. In the presented example it is not necessary, but might be of interest in case of noisy data.

• This is a nice approach -thank you. I have added your plots (which will only appear when peer reviewed). I regret I don't agree with your analysis. The plots for model 1 look nice and quadratic (as I suggested in the second part of my question). The plots for model 2 are wild and could easily take you off in the wrong direction. I would suggest that the comments you made can be reversed. The issue of parameters being of the same order of magnitude is one I will think about. – Hugh Aug 28 '15 at 17:07
• A very clear demonstration of going in circles. I wonder if this is almost a bug because it should be detectable and I would imagine there is a work around. The point about the phase is worth further thought. The result should be periodic in phase so there are parallel good valleys and it will not matter which one you are in. Thanks – Hugh Aug 31 '15 at 15:43
• @Hugh At other occasions Mathematica detects oscillatory behaviour, so it should be possible here as well. Overall, I think this is a good example that making a proper transformation of parameter space can remove pathological behaviour and result in much better conversion. The idea is to decouple the parameters and make things orthogonal such that you get a proper parabolic minimum. – mikuszefski Aug 31 '15 at 16:30
• I agree about trying to find a good parameter space. The counterintuitive part is that it is better to have two nonlinear parameters and one linear rather than two linear and one nonlinear. I had thought that the more linear parameters the better but that is wrong. I wonder if the final polishing is in a quadratic well. – Hugh Aug 31 '15 at 16:36
• What? No buttocks? – Dr. belisarius Sep 1 '15 at 17:03

It a small help, but I also did some analyses on this interesting problem and I found that the methods "NMinimize" and "ConjugateGradient" of NonlinearModelFit can find the real solution when the other methods fail, at least on MMA 10.2 (while on MMA 9 this did not seem to me to be effective).

With respect to the difficulty of getting the frequency right, the chart below shows the model error as a function of amplitude and f (the axis in front, between 0 and 1).

.

The fitting value should be 1/8 and it can be seen that the error decreases only for f values in a very thin strip around the optimal value. A little away, the error is almost constant. Consequently, only a thorough examination of the error behavior can point to the correct frequency, otherwise it will not be found, and amplitude and f will be practically arbitrary.

My suggestion is to fix the correct frequency with a Fourier analysis before any attempt to fitting the other parameters.

EDIT

The code to obtain the 3d plot was:

m1 = A Sin[2 \[Pi] f t] + B Cos[2 \[Pi] f t];
m2 = A Cos[2 \[Pi] f t + B];

data = Block[{A = 43, B = 94, f = 1/8},
Table[{t, m1 + 5 RandomReal[II]}, {t, 0, 30, 0.1}]
];
gdata = ListLinePlot[data]


fit1a = NonlinearModelFit[data, m1, {A, B, f}, t, Method -> Automatic];
Show[gdata, Plot[fit1a[t], {t, 0, 30}, PlotStyle -> Red]]


fit1b = NonlinearModelFit[data, m1, {A, B, f}, t,
Method -> "NMinimize"
];
Show[gdata, Plot[fit1b[t], {t, 0, 30}, PlotStyle -> Red]]


{x,y} = Transpose[data];
My = Mean[y];
ClearAll[error];
error[model_, {AA_, BB_, ff_}] :=
Block[{A = AA, B = BB, f = ff},
Norm[y - model /. t -> x]/My
];


The error plot above comes from:

Plot3D[error[m2, {103.3, B, f}], {B, -\[Pi], \[Pi]}, {f, 0, 1},
PlotRange -> All, AxesLabel -> {"B", "f"},
ColorFunction -> "Rainbow", MaxRecursion -> 3
]


With respect to Fourier:

afy = Abs@Fourier[y];
Pick[Range[0, Length[afy] - 1], afy, Max[afy]] // N


{4,297}

Note that 30 is the sampling time of the signal:

m1f = A Sin[2 \[Pi] (4/30 + f) t] + B Cos[2 \[Pi] 4/30 + f) t]

fit1f = NonlinearModelFit[data, m1f, {A, B, {f, 0}}, t,
Method -> Automatic
];
Show[gdata, Plot[fit1f[t], {t, 0, MT}, PlotStyle -> Red]]


Now, even with the Automatic method, since the frequency is centered around a close-to-optimal value, the fit comes out correctly.

• Thanks for looking at the problem. Please could you post your code so that we can fully see what you are doing? Your suggestion of finding the frequency from Fourier is extremely difficult as this post shows. The method here is an attempt to avoid using Fourier. – Hugh Sep 2 '15 at 21:18
• @Hugh, I tried to clarify above. The final step fits with the parameter f starting from zero since, in the frequency-adapted model m1f, the parameter f is centered around the value coming from the Fourier analysis. – user8074 Sep 2 '15 at 21:54

In short the two models presented are identical for the particular dataset.

"Intuition is sometimes tricky on fitting procedures." in the answer by @mikuszefsky is right on the money. So is the comment by @MarcoB about increasing the number of iterations.

The two fits are found just fine by NonlinearModelFit:

data = With[{a = 43.45582489316203,
f = 94.92003941300389, ϕ = 431.155471523826},
SeedRandom[1234];
Table[{t, a Cos[2 π f t + ϕ] + RandomReal[{-0.1, 0.1}]},
{t, 13.439999656460714, 13.479799655455281, 0.0002}]];
m1 = a Cos[2 π f t] + b Sin[2 π f t];
m2 = a Cos[2 π f t + ϕ];
nlm1 = NonlinearModelFit[data, m1, {{a, -30}, {b, 30}, {f, 95}}, t, MaxIterations -> 1000]
nlm2 = NonlinearModelFit[data, m2, {{a, 40}, {f, 100}, {ϕ, 3.9}}, t]


Both models are essentially identical but one cannot tell that from looking at the ParameterConfidenceIntervalTable.

The estimated variances are the same:

nlm1["EstimatedVariance"]
(* 0.0037219 *)
nlm2["EstimatedVariance"]
(* 0.0037219 *)


The $$AIC_c$$ values are identical:

nlm1["AICc"]
(* -543.167 *)
nlm2["AICc"]
(* -543.167 *)


The predictions are identical:

ListPlot[{Transpose[{nlm1["PredictedResponse"],
nlm2["PredictedResponse"]}], {{-45, -45}, {45, 45}}}, Joined -> {False, True}]


For m1 the numerical instabilities can be also seen by examining the parameter correlation matrix:

nlm1["CorrelationMatrix"]//MatrixForm


When the values in the off-diagonal are this close to 1 or -1, that suggests that there will either be numerical instability in the estimates and/or over-parameterization of the model for the particular data set.

The correlation matrix for m2 is also not without issues:

nlm2["CorrelationMatrix"]//MatrixForm


Given the values of the parameters for one model, one can find the parameters for the other model that give identical predictions. Consider the following slight deviation from the original definitions:

m1 = a Cos[2 π f t] + b Sin[2 π f t];
m2 = c Cos[2 π f t + ϕ];


(Only m2 is changed: a is replaced with c.) Then

{a, b} = c {Cos[ϕ], -Sin[ϕ]}


and

• Thanks for looking at this. I know the models are identical. My issue is why is a model with two nonlinear parameters better than one with one nonlinear parameter? How should we educate our intuition if this happens? Does make life difficult. – Hugh Aug 24 '20 at 20:04
• The simple reason is that the estimators for the parameters in m2 are less correlated than the estimators of the parameters in m1 (and that all of the correlations are so extreme in m1). That wrecks havoc with standard methods of optimization. The optimal solution either takes longer (more iterations) or a local optimum gets picked rather than the global optimum when the correlation is high. And it also depends on the available data. It's essentially the same thing about fitting a polynomial straightaway or using orthogonal polynomials. Same predictions but one more stable. – JimB Aug 24 '20 at 20:46
• This is the first comment that makes simple sense to me. I had not come across the idea that fitting is easier with less well correlated parameters. Makes good sense. Also, comment on orthogonal polynomials now makes sense. Generalising; transforming the data or the model to be as orthogonal as possible must be a good thing. Thank you. – Hugh Aug 25 '20 at 7:49