How does Mathematica compute the covariance matrix using "MultiNonlinearModelFit"?

Question

I am doing a simple multi-non linear fit. I get the covariance matrix. Does anyone know how Mathematica calculate this Matrix? I have found the computation when mathematica fit only one signal but have not found any definition when I perform the multi signals case.

ClearAll["Global`*"]
time = {0, 1, 3, 5, 6, 7};
For[k = 1, k <= 10, k++,
data = Flatten[Table[{Exp[1.8 1 + 1.1   x] + 
  RandomVariate[NormalDistribution[0, 10]]}, {x, 1, 6}]];
data2 = Flatten[Table[{Exp[2 1 + 1.1 1.2 x] + 
  RandomVariate[NormalDistribution[0, 10]]}, {x, 1, 6}]];
]
model[a, b] = Exp[1.8 a + b x^2];
model2[a, b] = Exp[2 a + 1.2 b x^2];
Listdata = Table[{time[[i]], data[[i]]}, {i, 1, Length[time]}];
Listdata2 = Table[{time[[i]], data2[[i]]}, {i, 1, Length[time]}];
dataa = Join[{Listdata}, {Listdata2}];
fit = ResourceFunction["MultiNonlinearModelFit"][dataa, {model[a, b], model2[a, b]}, {a, b}, {x}, Method -> "LevenbergMarquardt", MaxIterations -> Infinity];
param = fit["BestFitParameters"]
covmat = fit["CovarianceMatrix"] // MatrixForm

Thanks, Edoardo

Answer 1

The short answer is that the covariance matrix is estimated as if there was a single dataset with dummy variables indicating which subset of data belongs to which model.

MultiNonlinearModelFit is excellent and very convenient and the only potential issue I have with it is that it makes the unstated assumption that there is a common error variance which I would argue is rare to be true in practice.

So here is essentially what MultiNonlinearModelFit does:

Generate a data set where each subset of data is given a dummy variable that identifies the subset of data. For your example, I just use the values 0 and 1.

d0 = Transpose[Join[{ConstantArray[0, Length[Listdata]]}, Transpose[Listdata]]];
d1 = Transpose[Join[{ConstantArray[1, Length[Listdata]]}, Transpose[Listdata2]]];
(d = Join[d0, d1]) // MatrixForm

Now run NonlinearModelFit:

nlm = NonlinearModelFit[d, (1 - c) model[a, b] + c model2[a, b], {a, b}, {c, x}];
nlm["CovarianceMatrix"] // MatrixForm

This is exactly what you get from MultiNonlinearModelFit. If you want to know how NonlinearModelFit estimates the covariance matrix, that is a different question.

I repeat that this function should only be run if one knows that the error variance is the same among the models being considered. If one doesn't know, then at minimum taking a look at the residuals for each dataset is essential. (And one should look at the residuals in any event.)