Dear community, Does anyone know how Mathematica compute the unbiased variance (and the covariance matrix) in the case of a simultaneous fit of two signals (i.e. command "MultiNonlinearModelFit ")?

Hereby there is a working example in the case of a fit of a single signal:

ClearAll["Global`*"]


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"], \[Sigma] -> nlm["EstimatedVariance"]^0.5]
(*{a->1.5063204891556876`,b->1.4263297975161129`,\[Sigma]->0.\
810936632925338`}*)
nlm["CovarianceMatrix"]
(*{{1.2134907415092202`,-0.12830437622650265`},{-0.12830437622650268`,\
0.3601478400755359`}}*)



(*Construct the log likelihood function*)
logL = LogLikelihood[NormalDistribution[0, \[Sigma]], 
   data[[All, 2]] - Log[a + b data[[All, 1]]^2]];

(*Get maximum likelihood estimates*)
sol = NMaximize[{logL, \[Sigma] > 0 && a > 0 && b > 0}, {a, 
   b, \[Sigma]}]
(*{-6.039843700423039`,{a->1.5063204889322932`,b->1.4263297932385333`,\
\[Sigma]->0.6621269804210389`}}*)


(*Adjust estimator of \[Sigma] so that the estimator of \[Sigma]^2 is \
unbiased*)
n = Length[data];
p = 2;(*Number of fixed parameters to be estimated:a and b*)
sigma = \[Sigma] Sqrt[n/(n - p)] /. sol[[2]];
(*0.8109366234806664`*)
(*Modify the maximum likelihood estimate for \[Sigma]*)
sol[[2, 3]] = \[Sigma] -> sigma

Hereby there is an example in the case of a fit of two signal:

ClearAll["Global`*"]
time = {0, 1, 3, 5, 6, 7};
Clear[k]
For[k = 1, k <= 10, k++,
 data = Flatten[
   Table[{Log[1.8 1 + 1.1   x^2] + 
      RandomVariate[NormalDistribution[0, 10]]}, {x, 1, 6}]];

 data2 = Flatten[
   Table[{Log[2 1 + 1.1 1.2 x^2] + 
      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 = TimeConstrained[
   ResourceFunction["MultiNonlinearModelFit"][
    dataa, {model[a, b], model2[a, b]}, {a, b}, {x}, 
    Method -> "LevenbergMarquardt"], 10];
mle = Append[
  fit["BestFitParameters"], \[Sigma] -> fit["EstimatedVariance"]^0.5]

(*Construct the log likelihood function*)
logL = LogLikelihood[NormalDistribution[0, \[Sigma]], 
   data[[All]] - Log[1.8 a + b time[[All]]^2] + data2[[All]] - 
    Log[2 a + 1.2 b time[[All]]^2]];

(*Get maximum likelihood estimates*)
sol = NMaximize[{logL, \[Sigma] > 0 && a > 0 && b > 0}, {a, 
   b, \[Sigma]}] 

(*Adjust estimator of \[Sigma] so that the estimator of \[Sigma]^2 is \
unbiased*)
n = Length[data];
p = 2;(* Number of fixed parameters to be estimated:a and b *)
sigma = \[Sigma]  Sqrt[n/(n - p)] /. sol[[2]]

(*Modify the maximum likelihood estimate for \[Sigma]*)
sol[[2, 3]] = \[Sigma] -> sigma;

Thanks in advance