FindFit works, NonlinearModelFit doesn't

Question

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

Answer 1

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.