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

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

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

• Agreed that the unstated assumption about common variance can be dangerous. I should probably make a note about that the next time I work on it. That said: people also routinely throw NonlinearModelFit at data where this assumption is very obviously violated. Oct 18 '21 at 20:36
• @SjoerdSmit Yes, violating, ignoring, ignorance of assumptions in statistical analyses is rampant. I do wish Mathematica would include the capability for multiple sources of error (i.e., mixed models) in the high level functions. But I also understand that the programming of even linear mixed models is difficult enough due to convergence issues (many times caused by wrong models or bad data or both). And to repeat, your MultiNonlinearModelFit` is really well put together. Not using it and doing that stuff "by hand" is fraught with danger.
– JimB
Oct 18 '21 at 21:26
• Dear @JimB , Thank you very much. In the case I have 3 signals instead of 2, what is the dummy variable and how does the NonlinearModelFit change? Finally, do you know if there is something in the literature regarding this implementation? Thanks, Edoardo Oct 19 '21 at 7:53