How does Mathematica calculate the Standard error during a non linear fit?

Question

I am doing a simple non linear fit. If I use

nlm["ParameterTable"]

I get a table with Standard Errors. Does anyone know how Mathematica calculate those "standard errors"? I have searched but have not found any definition in the documentation.

Thanks, Umberto

Answer 1

Rather than give formulas here is how to duplicate the results (coefficient estimates and standard errors) from NonlinearModelFit.

First show the example from NonlinearModelFit:

data = {{0, 1}, {1, 0}, {3, 2}, {5, 4}, {6, 4}, {7, 5}};
nlm = NonlinearModelFit[data, Log[a + b x^2], {a, b}, x];
mle = Append[nlm["BestFitParameters"], σ -> nlm["EstimatedVariance"]^0.5]
(* {a -> 1.5063204891556876`,b -> 1.4263297975161129`,σ -> 0.810936632925338`} *)
nlm["CovarianceMatrix"]
(* {{1.2134907415092202`,-0.12830437622650265`},{-0.12830437622650268`,0.3601478400755359`}} *)
nlm["ParameterErrors"]
(* {1.1015855579614415`,0.6001231874169968`} *)

Now get the estimates from the use of the LogLikelihood function:

(* Construct the log likelihood function *)
logL = LogLikelihood[NormalDistribution[0, σ], 
   data[[All, 2]] - Log[a + b data[[All, 1]]^2]];

(* Get maximum likelihood estimates *)
sol = NMaximize[{logL, σ > 0 && a > 0 && b > 0}, {a, b, σ}]
(* {-6.039843700423039`,
   {a -> 1.5063204889322932`,b -> 1.4263297932385333`,σ -> 0.6621269804210389`}} *)


(* Adjust estimator of σ so that the estimator of σ^2 is unbiased *)
n = Length[data];
p = 2;  (* Number of fixed parameters to be estimated: a and b *)
sigma = σ Sqrt[n/(n - p)] /. sol[[2]]
(* 0.8109366234806664` *)
(* Modify the maximum likelihood estimate for σ *)
sol[[2, 3]] = σ -> sigma;

(* Now get the standard errors for the estimators of a and b *)
(* First construct the second derivative of logL symbolically *)
symbolicData = Table[{x[i], y[i]}, {i, n}];
logL = LogLikelihood[NormalDistribution[0, σ],
   symbolicData[[All, 2]] - Log[a + b symbolicData[[All, 1]]^2]];
h = D[logL, {{a, b}, 2}];
(* Now determine the expectation of h *)
Eh = h /. y[i_] -> Log[a + b x[i]^2] /. Table[x[i] -> data[[i, 1]], {i, n}];

(* Take the minus the inverse and plug in the maximum likelihood estimators to
   get the estimates of the variances and covariance of the parameter estimators *)
cov = -Inverse[Eh] /. sol[[2]]

(* Standard error for the estimator of a *)
cov[[1, 1]]^0.5
(* 1.1015855579614415` *)
(* Standard error for the estimtor of b *)
cov[[2, 2]]^0.5
(* 0.6001231874169969` *)