Obtain reduced chi-squared for fit

Question

Dear Colleagues, I have performed a fit according to the code below. The last step is to obtain the value for the reduced chi-squared for the goodness of fit. Does mathematica have a simple function that will give the answer directly or do I have to write the code to manipulate the data myself?

ClearAll["Global`*"]
E0 = 6   ;                                               
(* Electron energy in MeV*)

r0 = 4.975 ;   (* 10 MeV csda range in water cm*)
mc2 = 0.51099895 ; (* Electron mass*c2 MeV*)
q = 3 ; (* Electron Charge compatible with units of 
Mev, cm and teslas *)
lambda = E0/mc2 ;               (*constant*)
m0 = 0.0186269157429903;                                          

(*All*)
Zf[csda_] := (-1 + Exp[m0*csda])/m0;
arc[B_, mg_, rho_] := 
mg*(rho - (2*E0/(q*B))*ArcSin[rho/(2*E0/(q*B))]); 
NZfprime[mg_, rho_, B_] := Exp[m0*rho + arc[B, mg, 
rho]];
alpha[mg_?NumericQ, csda_?NumericQ, B_?NumericQ] := 
NIntegrate[NZfprime[mg, rho, B], {rho, 0, 
csda}]/Zf[csda];
data = {{0.3, 0.999734507157122}, {0.4, 
0.999385602795506}, {0.5, 
0.999000073501766}, {0.6, 0.99853529076513}, {0.7, 
0.998014296669706}, {0.8, 0.997441248453792}, {0.9, 
0.996780680891339}};
error = {1.48598975699963*10^(-4), 1.51681330269219*10^(-4), 
1.50344266461607*10^(-4), 1.58948671957842*10(^-4), 
1.56780798669199*10^(-4), 1.61928577268932*10^(-4), 
1.69830239808624*10^(-4)};

limit = Exp[rho*m0];
alpha[mg_?NumericQ, csda_?NumericQ, B_?NumericQ] := 
NIntegrate[If[B == 0, limit, NZfprime[mg, rho, B]], 
{rho, 0, csda}, 
Method -> {"GlobalAdaptive", Method -> 
"MultipanelRule"}]/Zf[csda];
solution = 
NonlinearModelFit[
data, {alpha[mg, csda, B] , mg > 0, 
0 < csda}, {{mg, 0.40}, {csda, 1.}}, B, Weights -> 
1/error^2, 
Method -> "NMinimize"]
(*{csda->0.162831,l->0.934341}*)
Show[Plot[solution[x], {x, 0, 1}], ListPlot[data]]

I understand that the question may seem trivial but I have been unable to find the function myself.

Regards, Ricardo

Answer 1

(* Reduced chisquare statistic:  the "hard way" *)
nObs = Length[data]
nParms = 2   (* mg and csda *)
r = solution["FitResiduals"];
r.Inverse[DiagonalMatrix[error^2]].r/(nObs - nParms)
(* 0.09953665214762722 *)

(* Reduced chisquare statistic: the "easy way" *)
solution["EstimatedVariance"]
(* 0.09953665214762722 *)

Just note that different fields use different names for this summary statistic. See Reduced chi-squared statistic. (The need to make up different names for technical terms from the originating field must be an interesting study in sociology and/or psychology. Consider that folks can't even be consistent on chi-square, chisquare, or chi-squared.)