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

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

## 1 Answer

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

• Thanks a lot. Could not hope for a better answer! Jul 7, 2017 at 8:49
• Yes, as @Umberto said, one could not hope for a better answer! Apr 13, 2018 at 14:41