# FindFit works, NonlinearModelFit doesn't

I have a problem using NonlinearModelFit vs FindFit. Here is my code.

vsGdata={{405, 1252.43},{404.022, 1593.51},{401.078,1027.},{396.132,1395.1},{389.123,1432.3},{379.959,1680.32},{368.512,1953.14},{354.609,2294.16},{338.013,2870.8},{318.399,4929.35},{295.312,65586.6}}
h = 6.62*10^-34;
c = 3*10^8;
e = 1.6*10^-19;
k = 1.38*10^-23;
PlanckPhotonDensity2 = a (2 e^3 10^-12 En^2)/(h^2 c^2) 1/(Exp[( 10^-6 e En)/(k T1)] - 1) + b (2 e^3 10^-12 En^2)/(h^2 c^2) 1/(Exp[( 10^-6 e En)/(k T2)] - 1);
func2[T1_?NumberQ, a_?NumberQ, T2_?NumberQ, b_?NumberQ, Enl_?NumberQ] := NIntegrate[PlanckPhotonDensity2, {En, Enl, \[Infinity]}]
fit = FindFit[vsGdata, func2[T1, a, T2, b, Enl], {{T1, 0.6}, {a, 1.3*10^18}, {T2, 0.06}, {b, 1.3*10^42}}, Enl]
nlm = NonlinearModelFit[vsGdata, func2[T1, a, T2, b, Enl], {{T1, 0.6}, {a, 1.3*10^18}, {T2, 0.06}, {b, 1.3*10^42}}, Enl]


Find fit returns

{T1 -> 0.643718, a -> 7.13416*10^17, T2 -> 0.0631721, b -> 2.04144*10^41}


which is a good fit.

NonlinearModelFit returns

FittedModel[func2[0.643718, 7.13416*10^17, 0.0631721, 2.04144*10^41]]


and then nlm["ParameterTable"] returns repeated

integrand ... has evaluated to non-numerical values for all sampling points in the region with boundaries ...


messages, and then

Missing[]


Why is there this difference between FindFit and NonlinearModelFit here, and is there a way to fix it? I'd really like to know the error on the fit...

• Running your code as is I get "Failed to converge to the requested accuracy or precision within 100 iterations" for both FindFit and NonlinearModelFit (Mathematica 10.4.1, Windows 10). I suspect you'll need to do some scaling of your data to get convergence. – JimB Jun 14 '16 at 17:28
• @Jim I don't think that message matters very much. I get a good fit from FindFit, and NonlinearModelFit seems to find the same parameters, but not return all properties. – NLambert Jun 15 '16 at 8:15

Your function is the sum of two functions of a similar form. The single function that covers both is

f[En_?NumberQ, a_?NumberQ, T_?NumberQ] :=
NIntegrate[(a e^3 x^2)/(500000000000 c^2 (-1 + E^((e x)/(1000000 k T))) h^2),
{x, En, ∞}]


I also scaled the parameters to be of similar size (and used essentially your estimated values for the starting values):

nlm = NonlinearModelFit[vsGdata,
f[Enl, a 10^17, T1] + f[Enl, b 10^39, T2],
{{T1, 0.6}, {a, 7}, {T2, 0.06}, {b, 3}}, Enl, MaxIterations -> 250];

sol = nlm["BestFitParameters"]
(* {T1 -> 0.6972919307351225,a -> 3.800671958879103,
T2 -> 0.0685099012593085,b -> 2.775150812998223} *)


Looking at the estimated correlations among the estimators one finds very, very high correlations:

nlm["CorrelationMatrix"] // MatrixForm


Such high correlations many times suggest an over or poorly parameterized model (which sometimes can be fixed with lots more data). (This is not to say that the 11 data points were not expensive or easy to obtain.) And such high correlations can also result in numerical issues (lack of convergence, etc.).

So where is the problem with how you originally did things? Here is a plot of the two functions along with their sum and the data using the estimated coefficients:

Show[ListPlot[vsGdata, PlotRange -> All, AxesOrigin -> {290, 0}],
Plot[{f[En, a 10^17, T1] /. sol,
f[En, b 10^39, T2] /. sol,
(f[En, a 10^17, T1] /. sol) + (f[En, b 10^39, T2] /. sol)},
{En, 295, 405}, AxesOrigin -> {290, 0}, PlotRange -> All],
ImageSize -> Large]


It looks like one of the functions only influences the sum for the first two data points and the remaining data points look pretty linear. (Again, this suggests over-parameterization given just 11 data points.)

While it might be a bit better to attempt to fit the log of the sum of the functions, I don't think it will change things much.