2
$\begingroup$

I want to understand the rule under which the function Nonlinearfit stops to give the final result. I know that if the number of iterations exceeds MaxIterations it stops, but otherwise it should give the best fit of the data for the given model.

In my experience instead it stops sometimes far before finding the best fit. One proof is that increasing the number of fitting parameters the quality of the fit sometimes decreases (often dramatically).

Consider for example the data in http://pastebin.com/raw.php?i=KqRHKE0pqþ and the program

minN = 1; maxN = 8;

g[x_,xo_,σ_,a_]:=a Exp[-((x-xo)^2/(2 σ^2))]

kvar[k_Integer]:=ToExpression@Map[StringJoin[#,ToString[k]]&,{"x","σ","a"}]

gmodel[n_Integer] := const + Sum[g[x, Sequence @@ kvar[i]], {i, 1, n}]

gpars[n_Integer]:=Flatten@{const,Array[kvar,n]}

myFittedModelsU1={};
Do[{
Print[kT];
tmpM=NonlinearModelFit[myData,gmodel[kT],gpars[kT],x,Method->{NMinimize,Method->{"DifferentialEvolution","ScalingFactor"->0.9,"CrossProbability"->0.7,"PostProcess"->{FindMinimum,Method->"QuasiNewton"}}} ];
AppendTo[myFittedModelsU1,tmpM];
},{kT,minN,maxN,1}]

The program tries to fit the data with model given by a mixture of Gaussian whose number changes for minN to maxN.

Increasing the number of parameters should increase the agreement between the model and the data (that does not mean the quality of the fit, of course), but for example in this case the fit is better for N=7 gaussians than for 10. Why this happens and how could I avoid this effect and obtain the best fit for my model?

$\endgroup$
  • 2
    $\begingroup$ I cannot run the code, but I do not see that you provide initial conditions. According to the documentation that should be done. In fact, I myself also fit many data that appears to be pretty simple but still fails if the initial values are not provided. You should definitely provide them if their values will be quite away from 1. $\endgroup$ – Lukas Oct 6 '15 at 19:09
  • $\begingroup$ Your pastebin link doesn't work. $\endgroup$ – MarcoB Oct 6 '15 at 19:11
  • 1
    $\begingroup$ I've been able to download the data and can certainly verify that the process stops too soon. But looking at the coefficients even for the 7 Gaussian terms, some of the multiplicative coefficients are negative (meaning that you're adding in upside down Gaussians) and others seem "extreme". While there are likely some techniques (such as standardizing both the dependent and independent variable) that would reduce the apparent numerical instability, I wonder if a small set of Gaussian curves (say < 10) can provide an adequate fit. What about nonparametric regression such as loess? $\endgroup$ – JimB Oct 7 '15 at 6:57
  • 1
    $\begingroup$ You have some fitparameters that you pass to NonlinearModelFit. Together with them you give initial conditions. It is right in the doc and easy to see. I suggest that you have a look at the whole doc on fitting... Seema you didn't do that $\endgroup$ – Lukas Oct 7 '15 at 17:14
  • 1
    $\begingroup$ I only played around a little bit with your fit, but saw a significant improvement when increasing WorkingPrecision. Example for kT=3 $\endgroup$ – Karsten 7. Oct 8 '15 at 13:38
3
$\begingroup$

This time a real answer. While looking at the estimated coefficients for the models that should have been better than the models with fewer Gaussians, it appeared that several of the fitted Gaussians were essentially just adding a constant rather than a normal shaped curve. So (roughly) the intercept and several of the "flat" Gaussians were duplicating each other. By removing the intercept and only allowing the multiplicative coefficients (a1, a2, ...) to be positive by putting those symbols into Exp[] (rather than putting in explicit constraints), the models seem to converge and fit properly and didn't need any further tweaking other than raising the maximum number of iterations.

To select the "optimal" number of Gaussians, I suggest choosing the number of Gaussians that minimized the AIC statistic. In the code below the number of Gaussians, mean square error, and AIC statistic are printed. For this data it appears that fitting more than 8 Gaussians are warranted.

(* Minimum and maximum number of Gaussians to consider *)
minN = 1; maxN = 8;

(* Define some functions to make life easier *)
g[x_, xo_, σ_, a_] :=  Exp[a - ((x - xo)^2/(2 σ^2))]
kvar[k_Integer] := ToExpression@Map[StringJoin[#, ToString[k]] &, {"x", "σ", "a"}]
gmodel[n_Integer] := Sum[g[x, Sequence @@ kvar[i]], {i, 1, n}]
gpars[n_Integer] := Flatten@{Array[kvar, n]}

(* Fit the models and print out the mean square error and AIC statistic *)
myFittedModelsU1 = {};
Do[{tmpM = NonlinearModelFit[myData, gmodel[kT], gpars[kT], x, MaxIterations -> 10000];
   AppendTo[myFittedModelsU1, tmpM];
   Print[{kT, myFittedModelsU1[[kT]]["EstimatedVariance"], myFittedModelsU1[[kT]]["AIC"]}]},
  {kT, minN, maxN, 1}]

(* Show the fits *)
Show[{ListPlot[Table[{x, myData[[x]]}, {x, Length[myData]}]],
  Plot[{myFittedModelsU1[[1]][y], myFittedModelsU1[[2]][y], 
        myFittedModelsU1[[3]][y], myFittedModelsU1[[4]][y],
        myFittedModelsU1[[5]][y], myFittedModelsU1[[6]][y],
        myFittedModelsU1[[7]][y], myFittedModelsU1[[8]][y]},
       {y, 0, Length[myData]},
       PlotStyle -> {Black, Red, Blue, Green, Orange, Gray, Cyan, {Purple, Thickness[0.005]}},
       PlotLegends -> {"1 Gaussian", "2 Gaussians", "3 Gaussians", 
       "4 Gaussians", "5 Gaussians", "6 Gaussians", "7 Gaussians", "8 Gaussians"}]}]

(* {1,4277.44,15011.6} *)
(* {2,2136.91,14084.7} *)
(* {3,1853.15,13896.7} *)
(* {4,806.086,12784.2} *)
(* {5,171.139,10710.6} *)
(* {6,164.663,10661.9} *)
(* {7,163.8,10657.8} *)
(* {8,132.644,10378.} *)

Gaussian fits

| improve this answer | |
$\endgroup$
  • $\begingroup$ I've noted that moving the parameter a in Exp[] makes the process much faster but what seems to really increase the quality of the fit is the lack of the const parameter to be fitted (besides the gaussian mixture). What I don't like of this solution is that the fit ends to depend on the absolute values of the data (that is quite arbitrary here), ie if I add a constant value to all the data the quality of the fit changes; eg, if I subtract 10^3 to all data it stops to work - because the model assumes all positive values. I ask myself if there is an optimal way to "scale" the data then. $\endgroup$ – psmith Oct 8 '15 at 13:20
  • $\begingroup$ Equivalently, I can define gmodel[n_Integer] := k + Sum[g[x, Sequence @@ kvar[i]], {i, 1, n}] and fix the value of k for all the fitting procedure. Your solution corresponds to k=0. Now some value of k works better than others, some doesn't work at all (eg if k>600). How to choose the best value for k? $\endgroup$ – psmith Oct 8 '15 at 13:29
  • $\begingroup$ In any case, besides the previous comment, I have tried extensively now this solution and I did not find a great improvement with respect to the procedure exposed in the question. $\endgroup$ – psmith Oct 8 '15 at 17:03
  • $\begingroup$ It does work with this data set. Is it some other data set that it doesn't work on? If so, can you share that additional data set? If you have negative values, just add the minimum value to all observations. If the fit is adequate (by either assessing the fit visually or by achieving a desired and specified mean square error), then everything's fine. Also, I know you seem wedded to Gaussian kernels but other forms of non-parametric estimation (such as loess) don't have these issues. $\endgroup$ – JimB Oct 8 '15 at 17:19
  • 1
    $\begingroup$ I tried subtracting 0, 100, 200, 500, and 600 from the observed values with 8 Gaussians and the mean square error ranged from 132.644 (with subtracting 0) to 133.016 (with subtracting 600). That's a pretty small difference. I would strive for reasonable fits (using mean square error and residual plots as guides) rather than going after the optimal number of Gaussians. $\endgroup$ – JimB Oct 8 '15 at 20:55
4
$\begingroup$

Disclaimer: This is not a direct answer to your question but rather provides an alternative that still uses curves with a Gaussian shape.

While I think fitting such data is better described with nonparametric regression models (gams, loess, kernel regression, etc.), if I remember correctly, you were asked/commanded to use Gaussian curves.

As you've found despite setting some of the tuning parameters, NonlinearModelFit does not always automatically produce the fit that minimizes the sum of squares for this data and suite of models.

The alternative proposed is to use a fixed set of Gaussians equally-spaced as a set of basis functions to describe the data. This approach does end up with many more Gaussians to describe the data but convergence to the appropriate estimates is much more likely.

Below is some code that I've modified from your original code.

(* Min and Max of the predictor variable for this dataset *)
xmin = 1;
xmax = Length[myData];

(* Functions to create the model and list of parameters *)
delta[k_] := (xmax - xmin)/(k - 1);
g[x_, i_, k_, σ_, a_] := a Exp[-((x - (xmin + (i - 1) delta[k]))^2/(2 σ^2))]
kvar[k_Integer] := ToExpression@Map[StringJoin[#, ToString[k]] &, {"a"}]
gmodel[n_Integer] := a0 + Sum[g[x, i, n, σ, Sequence @@ kvar[i]], {i, 1, n}]
gpars[n_Integer] := Flatten@{σ, a0, Array[kvar, n]}

(* Estimate coefficients for 8 and 25 Gaussian basis functions *)
n = 8;
m8 = NonlinearModelFit[myData, gmodel[n], gpars[n], x, MaxIterations -> 5000];
n = 25;
m25 = NonlinearModelFit[myData, gmodel[n], gpars[n], x, MaxIterations -> 5000];

(* Display the mean square error *)
m8["EstimatedVariance"]
(* 1736.37 *)
m25["EstimatedVariance"]
(* 143.565 *)

(* Show fits *)
Show[{ListPlot[Table[{x, myData[[x]]}, {x, Length[myData]}]],
  Plot[{m8[x], m25[x]}, {x, xmin, xmax}, PlotStyle -> {Green, Red}, 
   PlotLegends -> {"8 Gaussians", "25 Gaussians"}]}]

Gaussian basis function fits

| improve this answer | |
$\endgroup$
  • $\begingroup$ I appreciate this solution, but I need a fit with the "optimal" number of gaussians. In this case you usually use many more than the necessary to obtain a good fit. $\endgroup$ – psmith Oct 7 '15 at 16:25
  • 1
    $\begingroup$ Understood. I did try feeding the estimates from the previous fit as starting values in your code (with a = 0 for the "new" Gaussian) but that didn't seem to help. Another grasping at straws thing I tried was to restrict some of the parameters to be positive (such as the "a") values but that, too, didn't seem to help. $\endgroup$ – JimB Oct 7 '15 at 16:33
  • $\begingroup$ I also tried to restrict all the parameters to "reasonable" values but, as you say, it does help much: still the estimated variance is far to be a monotonic function of the number of gaussian but instead oscillates in a quite chaotic way... $\endgroup$ – psmith Oct 7 '15 at 16:39

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.