Why is FindFit failing to find a suitable trig model?

Question

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.

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

Answer 1

@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}]]

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]}
 ]

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.

Answer 2

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}*)

Answer 3

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] 
Show[ListPlot[tidedata],
     Plot[Evaluate[model03[b] /. fit03[b /. fit04] /. fit04], {x, 0, 25}]
    ]  

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