5
$\begingroup$

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.

enter image description here

Thanks, Umberto

$\endgroup$
14
$\begingroup$

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` *)
$\endgroup$
  • 1
    $\begingroup$ Thanks a lot. Could not hope for a better answer! $\endgroup$ – Umberto Jul 7 '17 at 8:49
  • $\begingroup$ Yes, as @Umberto said, one could not hope for a better answer! $\endgroup$ – Anton Antonov Apr 13 '18 at 14:41

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.