stop rule for Nonlinearfit and best fit

Question

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?

Answer 1

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

Answer 2

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