# Problem with NonlinearModelFit

I'm having trouble with a non-linear fit:

fit = NonlinearModelFit[data, y0 + A Sin[\[Pi] (x - xc)/w], {y0, xc, A, w}, x]

where data has about 15 thousand points and looks like this:

Mathematica gives the following adjusted parameters:

{y0 -> 30.4428, xc -> 1.54318, A -> -0.000528519, w -> 0.999975}

However, this is a terrible fit:

I wouldn't complain, except that there is an obviously better fit:

with the parameters:

{y0 -> 30.45775, w -> 752.71185, A -> 3.62443, xc -> 872.72035}

This best fit is returned by Origin without effort. What can I do to achieve the same with Mathematica? (I like Mathematica!)

-
I understand that the question is not just about how to fit this particular dataset, but I can't help wondering, why are you fitting a sine function to this obviously not sine data? If you are trying to measure the period / amplitude, there are better ways to do that. –  Szabolcs Feb 25 '12 at 10:31
@Szabolcs I'm not just trying to get the period. My end purpose is to represent the data as a superposition of sines with distinct frequencies (which do not have to be all multiples of a fundamental frequencies, so it's not a Fourier series). On the other hand, I would like to know what are those better ways you mention to measure the period. –  becko Feb 25 '12 at 15:45
It's possible to get a good estimate by interpolating the Fourier transform. Even if the signal is a sum of sines with periods that are not multiples of a base frequency, the right way to do this is a Fourier transform. The base wavelength in the transform will be the length of the whole dataset, and the components will be integers fractions of that, so you get better resolution. Interpolating it (as in the link above) will let you detect wavelengths that are not integer fractions of that. –  Szabolcs Feb 25 '12 at 16:12
This was obtained by using Eureqa, a very nice and free tool –  belisarius Apr 3 '12 at 12:50

I have mentioned this in a comment already, but this seems like a good opportunity to provide some related discussion in the form of a full-fledged answer.

In Mathematica 8, we can take advantage of NMinimize to fit this data, using the Method -> NMinimize option of NonlinearModelFit. (This should also have worked in Mathematica 7, but unfortunately NMinimize was not recognised as a valid Method setting until version 8 due to a bug.) In particular, Storn-Price differential evolution, available to NonlinearModelFit using the option

Method -> {NMinimize, Method -> "DifferentialEvolution"}


has a lot to offer in this case, especially if you know a bit about how differential evolution works. This algorithm, as implemented in Mathematica, is documented at tutorial/ConstrainedOptimizationGlobalNumerical#24713453.

From the documentation, we see that the scaling factor $s$ (called $F$ by Storn and Price in their publication on the method and usually elsewhere) acts as an amplification factor on the scale of the global search. Thus, a large value of $s$ encourages more expansive searching of the parameter space, while small values encourage more intense exploration around local minima. Classically, $s$ can take values between 0 and 2, although Mathematica doesn't enforce this restriction. In practice one finds that values larger than unity cause an extreme expansion of the parameter space under search, which may be counterproductive. A "large" value of $s$, then, is something close to 1, and this is what we need in the current case since we may suspect that the initial values chosen for the parameters might be rather far from the global optimum, and do not want to risk falling into some local minimum along the way.

The behaviour of differential evolution with respect to crossover probability, $\rho$ (which, as pointed out by Daniel Lichtblau, is equal to Storn and Price's $1 - CR$), is also very important. Noting that two of the parameters, w and xc, are strongly correlated, and knowing that in such cases vigorous mutation is usually the most effective strategy, we might also consider setting $CR \approx 1$, i.e. $\rho \approx 0$. While the default value of $\rho = 0.5$ does work for this example, if more sine functions are introduced into the model, reducing $\rho$ will be practically mandatory.

Plenty of discussion (indeed, an extensive literature) on tuning the differential evolution parameters, including the (usually) less critical population size parameter, $m$ (a.k.a. $NP$), can be found elsewhere, if necessary. However, it's worth noting that the "correct" values may differ between Mathematica's implementation and others, especially for small populations, due to slight differences in the way that the three existing random points are chosen to produce new trial search points.

So, writing down our conclusions from the above, we have:

data = Import["dat.csv"]; (* with thanks to @Szabolcs *)

fit = NonlinearModelFit[
data,
y0 + A Sin[Pi (x - xc)/w],
{y0, xc, A, w}, x,
Method -> {NMinimize,
Method -> {"DifferentialEvolution",
"ScalingFactor" -> 0.9, "CrossProbability" -> 0.1,
"PostProcess" -> {FindMinimum, Method -> "QuasiNewton"}
}
}
]


Where one should note the undocumented in this context, albeit rather obviously existent, Method option for FindMinimum as used by NMinimize as used by NonlinearModelFit (yes, that's right: we are setting a Method's Method's Method!). This serves to hone the parameter values produced by differential evolution given that the latter is, by design in this case, not as efficient for local optimization as other methods. Here "QuasiNewton" corresponds to the method of Broyden, Fletcher, Goldfarb, and Shanno (BFGS), but "LevenbergMarquardt" could also have been used.

This gives us:

Or, as a list of rules:

{y0 -> 30.4578, xc -> 120.008, A -> 3.62443, w -> -752.712}


This a result consistent (up to the sign of w and the value of the phase factor xc) with that given by Origin. Was it achieved without effort (if this is considered important)? While this is inherently a subjective question, in my opinion, the answer is yes. No manually chosen initial values in sight!

A plot of the resulting model also makes it clear that this is a reasonable outcome (although evidently one could do better with a more involved model):

-
I think the only difference between our implementation and the Storn and Price algorithm was that we managed to reverse the sense of the crossover parameter. –  Daniel Lichtblau Feb 25 '12 at 20:28
@DanielLichtblau also, the difference vectors are chosen without any constraint, but in the Storn-Price version they're chosen to be mutually different as well as different from the target vector. This may bias the results for small populations and/or require slightly different values of $F$ and $CR$ to the ones normally recommended. But yes, the reversal of the meaning of $\rho$ vis-a-vis $CR$ is the most significant difference. And well noticed, as I missed this myself and wondered why $\rho \approx 1$ wasn't working "properly". –  Oleksandr R. Feb 25 '12 at 20:53
(I'm supposed to notice that sort of thing. Some of it may have been my code at one time. I actually only realized the reversal when I went to write up a paper last year and wanted to check that polarity. It matters more for discrete optimization. Am hoping Price and Storn don't get mad at me...) –  Daniel Lichtblau Feb 25 '12 at 21:29
@DanielLichtblau From my understanding, differential evolution isn't generally considered appropriate for discrete problems, as while it can do the job with a suitable objective function, it's very inefficient in terms of function evaluations. BTW, the rand/1 mutation strategy is alright, but others are IMO worth having too, as is an elitist selection method. (I've played around quite a bit with a DE implementation of my own in Python.) –  Oleksandr R. Feb 25 '12 at 22:07

Some time ago I wrote a tutorial on (very) basic fitting in Mathematica. Topics covered (or more like talked about) are error bars, how to construct fitted models, adding fancy confidence bands, and getting the fit parameters. You can download the file here. Maybe some of that is helpful to you.

Also, some eyecandy from the nb so people actually look at this answer:

(Gaussian fit to test data generated in the nb, plus confidence bands.)

-
that is a nice tutorial –  acl Feb 25 '12 at 2:11
Have to agree with @acl, nice tutorial, +1. –  rcollyer Feb 25 '12 at 3:28

I agree with @rcollyer that starting values are often necessary for non-linear model fits. Here is a nifty way (though slow) to use Manipulate to dynamically change starting values in and obtain some information about the fit and a plot of the model against the data.

Manipulate[mod = NonlinearModelFit[Transpose[{xx, y}],
y0 + A Sin[\[Pi] (x - xc)/w],
{{y0, y00}, {xc, xc0}, {A, A0}, {w, w0}}, x];
Show[
Plot[mod[t], {t, 0, 12000}, PlotRange -> {{0, 12000}, {25, 39}}],
ListPlot[Transpose[{xx, y}], PlotStyle -> Directive[PointSize[0]]],
PlotLabel ->
Column[{Row[{"AIC -> ", mod["AIC"]}],
Row[{mod["BestFitParameters"]}]}]],
{y00, -100, 100}, {xc0, 1, 1000}, {A0, -10, 10}, {w0, 1, 1000}
]


This gives you something like the following...

You can change the starting values to see just how dependent this particular fit is to their choice. The AIC value can be used to judge relative fit for each model (smaller is better).

It seems to me that the particular model you are trying to fit doesn't actually fit these data all that well...

NOTE: I used the results from my own answer to this question to obtain the data from the image provided by becko.

-
I was thinking about doing exactly that to get at his data. :) –  rcollyer Feb 25 '12 at 3:24
I think it would be simpler to do the Manipulate without the NonlinearModelFit, just plotting the data and the function y0 + A Sin(Pi*(x-xc)/w) and manipulating the parameters to see visually what fits better. Then we use those parameters as starting value for the Fit. That way the manipulate would be much faster. –  becko Feb 25 '12 at 15:59
However, my goal is to use a sum of various sines to fit the data, so in the end I'm going to have a lot of parameters. Is there any other way of increasing the accuracy of the fit without specifying initial values? –  becko Feb 25 '12 at 16:06
Manipulating the parameters is certainly an option but especially when there are many it can be very difficult to find ones that work well simultaneously. Though it slows things down, leaving the model fitting function in does much of the grid searching work for you. –  Andy Ross Feb 25 '12 at 18:02
Keep in mind that NonlinearModelFit is a very general function. If Origin has some dedicated function for fitting sine functions it is likely using Fourier methods like those mentioned by @Szabolcs. A general solver that doesn't take advantage of properties of the functional form is going to need good starting values. –  Andy Ross Feb 25 '12 at 18:28

I don't have time for a complete answer, but non-linear fits are inherently more difficult. So, you often need to specify initial values for your coefficients to get a good fit. This has the form:

NonlinearModelFit[data, y0 + A Sin[\[Pi] (x - xc)/w],
{{y0, yval}, {xc, xcval}, {A, Aval}, {w,wval}}, x]

-
What I don't understand is how Origin is able to do it without initial values. Origin has a fit specially designed for sine functions, a "sine fit". Maybe it exploits somehow the knowledge that the fit function is a sine. –  becko Feb 25 '12 at 16:02
Isn't there a different approach that doesn't involve specifying initial values? My goal is to express the data as a linear combination of sines, so I'm going to have a lot of parameters and I can't guess initial values for all. –  becko Feb 25 '12 at 16:03
@becko my internet connection has been intermittent since your comments, so my apologies for taking so long to reply. With your stated goal of using a linear combination of sines, I'd suggest using Szabolcs suggestions. Also, Oleksandr's method doesn't involve additional parameters whatsoever. –  rcollyer Feb 26 '12 at 2:27

There may be a better way of approaching this problem. Since what you are trying to do is to approximate a data set with a sinusoid, why not use the Fourier transform? After all, this is what it is intended for! Start by massaging the data a bit (shifting to the origin and centering) so that it is easy to see what is happening, and taking the Fourier transform:

dat = data[[All, 2]] - Mean[data[[All, 2]]];
fft = Fourier[dat, FourierParameters -> {-1, 1}];
ListPlot[Abs[fft[[1 ;; 20]]], PlotRange -> All, Filling -> Axis]


Since we are looking for a low frequency, I have just plotted the first few terms. Clearly the eighth term is the largest, and this gives the frequency of the best-fit sinusoid. This can be done automatically picking the location of the largest peak:

ind = First@First@Position[Abs[fft], Max[Abs[fft]]];
Show[ListPlot[Pi Abs[fft[[ind]]] Sin[ 2 Pi (ind - 1)/Length[dat]*Range[Length[dat]]],
PlotRange -> All, PlotStyle -> Red], ListPlot[dat]]


Note that we haven't used any "initial values", rather, the max in the FFT has shown us the best values. Moreover, since the OP says in a comment that "My goal is to express the data as a linear combination of sines" this is exactly what the FFT does -- simply take more terms and you will have a better approximation. Alternatively, if it is really desirable to do this using an optimization method, the FFT values could be used to provide good starting values for the optimization.

-