I'm trying to get FindFit (or NonlinearModelFit) to find a best fit trig function.

The result I achieve is so obviously not the good fit I'd like as you can see below. I think this should be simple, but I don't seem to be able to have it work correctly.

not findfit

Hopefully someone will be able to put me on the right track.

  • 3
    $\begingroup$ How is tidedata defined? $\endgroup$ – Coolwater Feb 2 '18 at 19:59
  • 3
    $\begingroup$ You need to choose appropiate starting values for the paramters. E.g. {a, -Subtract @@ MinMax[tidedata[[All, 2]]]} and {d, Median[tidedata[[All, 2]]]} $\endgroup$ – Coolwater Feb 2 '18 at 20:01
  • 5
    $\begingroup$ Do you have copy and paste-able code? $\endgroup$ – ktm Feb 2 '18 at 20:52
  • $\begingroup$ I would be interesting to have an answer to this question with the minimal guesses on the starting values or the contraints. I think one could obtain something like 1) no guesses concerning a,c,d 2) b has one contraint : expected result beetween 0 and 10 periods in the Plot. $\endgroup$ – andre314 Feb 2 '18 at 22:12
  • 3
    $\begingroup$ Please submit copyable data. $\endgroup$ – David G. Stork Feb 2 '18 at 22:46

@Coolwater gave the answer: one needs good starting values.

The default starting value of 1 for all parameters doesn't always work. While it would be more convenient if one wouldn't have to produce starting values, having some idea as to what the plausible values might be is usually a good idea. Fortunately the OP plotted the results and saw that there was a problem.

I "digitized" the data to reproduce the results:

tidedata = {{0.0322997, 2.45303}, {2.03488, 3.77871}, {3.03618, 3.98747},
 {4.00517, 3.8309}, {6.04005, 2.54697}, {8.01034, 1.24217}, {8.85013, 1.09603},
 {10.0452, 1.13779}, {11.9832, 2.37996}, {14.0504, 3.75783}, {15.1163, 4.07098},
 {15.9884, 3.87265}, {18.0556, 2.65136}, {19.9935, 1.26305}, {21.3178, 0.782881},
 {22.0284, 1.01253}, {23.9987, 2.31733}};

model = a Sin[b x + c] + d
fit = FindFit[tidedata, model, {a, b, c, d}, x]
(* {a -> -0.417056, b -> 0.894376, c -> 1.72554, d -> 2.48805} *)

If better starting values are given, then the appropriate estimates are found. It looks like there are two full cycles over a length of 25 which means that b could be estimated by $4\pi/25$ which is close to 0.5.

fit = FindFit[tidedata, model, {a, {b, 0.5}, c, {d, Mean[tidedata[[All, 2]]]}}, x]
(* {a -> 1.53685, b -> 0.519602, c -> -0.0399239, d -> 2.48761} *)
Show[ListPlot[tidedata], Plot[model /. fit, {x, 0, 25}]]

Better fit to data

Why the need for better starting values? In this case it's because the likelihood function is extremely "bumpy" and I'll show what is meant by that. Suppose we know a and d and have some idea of the standard error of estimate. Then we can look at the likelihood function for b and c.

Below I show that likelihood function is bumpy and the result depends on the starting values used.

(* Get some starting and ending values *)
start1 = {1, 1};
start2 = {1, -3};
start3 = {0.8, -2};
end1 = {b, c} /. FindFit[tidedata, model, {a, {b, start1[[1]]}, {c, start1[[2]]}, d}, x];
end2 = {b, c} /. FindFit[tidedata, model, {a, {b, start2[[1]]}, {c, start2[[2]]}, d}, x];
end3 = {b, c} /. FindFit[tidedata, model, {a, {b, start3[[1]]}, {c, start3[[2]]}, d}, x];

Now show the starting and ending values on a contour plot of the likelihood function:

σ = 0.1;
logL = LogLikelihood[NormalDistribution[0, σ],
   tidedata[[All, 2]] - (model /. {x -> tidedata[[All, 1]], 
       a -> 1.5368524967407906, d -> 2.4876143581877814})];
ContourPlot[logL, {b, 0.3, 1.1}, {c, -π, π}, 
 PlotLegends -> Automatic,
 Frame -> True, FrameLabel -> (Style[#, Bold, 18] &) /@ {"b", "c"},
 Contours -> -{100, 200, 300, 400, 500, 1000, 1500, 1600, 1700, 1800, 
    1900, 2000, 2400, 2500, 2600, 2700, 2800, 3000, 3500},
 Epilog -> {PointSize[0.025], Thick,
   Line[{start1, end1}], Red, Point[end1], Blue, Point[start1],
   Line[{start2, end2}], Red, Point[end2], Blue, Point[start2],
   Line[{start3, end3}], Green, Point[end3], Blue, Point[start3]}

Likelihood contour plot

The blue dots represent the 3 different starting values. The two red dots are the resulting wrong answers. The green dot represents the correct answer.

We see that starting values for 1 end up in a trough and the starting values for 2 end up on a local peak. The global peak is found when the starting values are a bit closer to the global peak and there aren't any bumps in the way.

It is simply a fact of life that sometimes really good starting values are needed. This is especially true for fitting multiple cycles of a sine wave.

| improve this answer | |

If you choose the "right" method, FindFit works without starting values:

model = a Sin[b x + c] + d
fit = FindFit[tidedata, model, {a, b, c, d}, x, Method -> "NMinimize"]
(*{a -> -1.53685, b -> 0.519602, c -> -3.18152, d -> 2.48761}*)

enter image description here

| improve this answer | |
  • 2
    $\begingroup$ +1 NMinimize is definitely more stable for such a function and seems to work even setting all of the default values to 1,000. However, a minor quibble about your wording. I think it would be more accurate to say "with the default starting values" rather than "without starting values". Of course, "without user-supplied starting values" would be good, too. $\endgroup$ – JimB Feb 3 '18 at 18:32
  • $\begingroup$ @JimB: Thanks. So let's call it "without user-supplied starting values". For me the main advantage of NMinimize is the abilility to achieve an "objective" minimum by only(if necessary) constraining the parameter range. $\endgroup$ – Ulrich Neumann Feb 4 '18 at 14:09

Experimental Approach.

Feel free to criticize, or to find this stupid, but with explanation it would be more welcomed.

Here is an approach that works without starting values for 3 of the 4 parameters.

The idea is to use the robustness of linear algebra, because the problem is in fact linear in 3 parameters :

It can be rewritten as : model = aa Sin[b x] + cc Cos[b x] + d

So one can write a FindFit over b that call another Findfit of a linear model over aa,cc,d.

Here is the code :

model03[b_]:= aa Sin[b x] + cc Cos[b x] + d  

fit03[b_?NumericQ]:= FindFit[tidedata, model03[b], {aa, cc, d}, x]  

fit04 = FindFit[tidedata, model03[b] /. fit03[b], {{b,.5
}}, x] 
     Plot[Evaluate[model03[b] /. fit03[b /. fit04] /. fit04], {x, 0, 25}]

enter image description here

Normally, the first FindFit (I mean the one over parameters aa,cc,d) doesn't use any "starting values" (in particular "default starting values")

There's a remaining question : what method is used by the FindFit over b ? (method probably changed due to the effect of ?NumericQ in fit03[b_?NumericQ]).

| improve this answer | |
  • $\begingroup$ Thanks to everyone - that helps enormously. $\endgroup$ – Tony Halsey Feb 5 '18 at 5:05

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.