# How does mathematica calculate CovarianceMatrix in NonLinearModelFit

I have read "StatisticalModelAnalysis-Tutorial" and try to follow the definitions to calculate the ParameterError step-by-step. What I figured out is ParameterErrors in nlm["ParameterTable"] is calculated from diagonal entries of Covariance Matrix Sqrt[nlm["CovarianceMatrix"]] // Diagonal. From link above the definition of CovarianceMatrix is as follows:

"The matrix is sigma^2 (X^T W X)^(-1), where sigma^2 is the variance estimate, X is the design matrix for the linear approximation to the model, and W is the diagonal matrix of weights."

This is the Function, data and weights I putted in NonLinearModelFit:

f = Function[{a, b, x}, b + a*x];
data = {{6., 21.}, {8., 31.}, {9., 39.}};
weights = {1., 1., 1.};
nlm = NonlinearModelFit[data, f[a, b, x], {a, b}, x, Weights -> weights, VarianceEstimatorFunction -> (1 &)];


The "DesignMatrix" X I think is produced from y_i estimates:

para = a /. nlm["BestFitParameters"];
parb = b /. nlm["BestFitParameters"];
yestimates = f[para, parb, data[[All, 1]]]
designmatrix = DesignMatrix[data, f[para, parb, x], x];
designmatrix // MatrixForm


Now one point I do not understand. What is sigma^2 ? I thought it could be EstimatedVariance:

varianceestimate = nlm["EstimatedVariance"];


Which in my case gives 1 because I did VarianceEstimatorFunction-> (1&) in nlm. The "Matrix of weights" W I thought to be:

matrixofweights = DiagonalMatrix[weights];


And now calculate the Covariance Matrix from definiton above:

varianceestimate*Inverse[Transpose[designmatrix].matrixofweights.designmatrix]//MatrixForm


out:{{6.08061, -0.189471},{-0.189471, 0.00624628}}

And comparison with original CovarianceMarix don't give the same result:

nlm["CovarianceMatrix"] // MatrixForm


out: {{0.214286, -1.64286}, {-1.64286, 12.9286}}

Can somebody help?

Edit: By the way, at least there is no factor between these different results. so my mistake could not only due to sigma^2.

• Are you sure you're not confusing Covariance and Variance? en.wikipedia.org/wiki/Covariance en.wikipedia.org/wiki/Variance – Feyre Aug 14 '16 at 12:41
• Parameters standard error is a standard deviation which is calculated from Sqrt[nlm["CovarianceMatrix"]] // Diagonal in this case I'm sure... I'am not sure in expression "variance estimate" – Schrubber Aug 14 '16 at 12:53
• So, now I nearly figured out myself. Should I edit the entry above or just do a comment? Sorry, for asking then solve it ... but I tried over hours don't come to a result.. – Schrubber Aug 14 '16 at 14:12
• Nothing wrong with that! Feel free to post an answer down the page to your own question. This sort of question doesn't often get an answer from the community because it's more statistics oriented than Mathematica oriented. I just hope next time others can be of more use. – Feyre Aug 14 '16 at 14:20
• so, where should or can I ask this kind of questions? – Schrubber Aug 14 '16 at 15:05

Just a note: it has been called the design matrix for linear models long before even I was born. A statistics class would help with the terminology.

However, in a nonlinear setting it turns out to be a slightly more complicated item that still results in the usual design matrix when the model is linear.

I'm going to ignore the issue of weighting because it's all about what X is for a nonlinear model. Below is shown an example of a nonlinear model using NonlinearModelFit and a brute force maximum likelihood method to show the underlying calculations.

 (* Generate data from a nonlinear model *)
SeedRandom[12472475];
x = Flatten[Table[{0, 1, 2, 3, 4, 5, 6, 7, 8}, {i, 5}]];
n = Length[x];
aa = 1;
bb = 0.5;
σσ = 0.05;
y = aa Exp[-bb x] +
RandomVariate[NormalDistribution[0, σσ], n];

(* Fit with NonlinearModelFit *)
nlm = NonlinearModelFit[Transpose[{x, y}], a Exp[-b z], {a, b}, z];

(* Fit with brute force method *)
(* Get log of the likelihood *)
logL = LogLikelihood[NormalDistribution[0, σ], y - a Exp[-b x]];
(* Find maximum likelihood estimates of parameters *)
mle = FindMaximum[{logL, σ > 0}, {{a, 1}, {b, 0.5}, {σ, 0.05}}];

(* Now approximate and then estimate the covariance matrix *)
(* Note the definition of X, the design matrix *)
X = D[a Exp[-b x], {{a, b}}];
(* Multiplication by n/(n-2) is to use the unbiased estimate of σ^2 *)
(* n - 2 is really sample size minus number of parameters *)
cov = (n/(n - 2)) (σ^2 Inverse[Transpose[X]. X]) /. mle[[2]];

(* Parameter estimates *)
nlm["BestFitParameters"]
(* {a -> 0.9842981255137611,b -> 0.5023511174255275} *)
mle[[2]][[1 ;; 2]]
(* {a -> 0.9842981255018496,b -> 0.50235111740454} *)

(* Estimates of covariance matrix *)
nlm["CovarianceMatrix"]
(* {{0.00034959357568038475,0.00016592363407431294},
{0.00016592363407431294,0.0002923502562857745}} *)
cov
(* {{0.00034959357700581665,0.0001659236346998225},
{0.0001659236346998225,0.00029235025737224096}} *)

Rationalize[cov/nlm["CovarianceMatrix"], 0.00001]
(* Estimates of σ^2 *)
nlm["EstimatedVariance"]
(* 0.002014630678417603 *)
((n/(n - 2)) σ^2 /. mle[[2, 3]])
(* 0.0020146306860807216 *)

• Thanks a lot. That help me to understand! – Schrubber Aug 16 '16 at 6:05

I will answer myself half the way. So I misunderstood the design matrix, it should be DesignMatrix OF X. ok that's stupid because it is called X. design matrix has the form {{1.,Xestimate_1} ,{1.,Xestimate_2},...}. The following code will produce the same matrix as nlm["CovarianceMatrix"] does. BUT just for linear models b + a*x. So there is still a open question. :-)

f = Function[{a, b, x}, b + a*x];
data = {{6., 21.}, {8., 31.}, {9., 39.}};
weights = {0.5, 4.5, 0.1};
nlm = NonlinearModelFit[data, f[a, b, x], {a, b}, x,Weights -> weights, VarianceEstimatorFunction -> Automatic]
yestimates = Table[nlm[data[[i, 1]]], {i, Length[data]}]
xestimates = Flatten[Table[x /. Solve[nlm[x] == yestimates[[i]]], {i,Length[data]}]]
xdesignmatrix = Table[{1., xestimates[[i]]}, {i, Length[data]}];
varianceestimate = nlm["EstimatedVariance"];
matrixofweights = DiagonalMatrix[weights];
varianceestimate*Inverse[Transpose[xdesignmatrix].matrixofweights.xdesignmatrix]
% // MatrixForm
nlm["CovarianceMatrix"] // MatrixForm


output: {{25.9587, -3.29752}, {-3.29752, 0.421488}}