Unexpected result computing covariance of complex data matrices?

Questions
1

It seems to me that Mathematica when computing the (sample) covariance for a complex valued data matrix computes it the wrong way!
If we have something like

Data = {{1.2 + 0.5 I, 2.6 - 3.1 I}, {5.3 - 0.5 I,1.1 + 5.6 I}, {3.5 + 3.5 I, 1.6 - 3.1 I}, {2.5 + 3.5 I, 4.6 - 3.1 I}, {4.3 - 0.5 I, 2.1 + 5.6 I}, {3.5 + 3.5 I, 4.6 - 3.1 I}, {4.3 - 0.5 I, 2.1 + 5.6 I}, {1.5 + 3.5 I, 3.6 - 3.1 I}, {2.3 - 0.5 I, 1.1 + 5.6 I}, {1.5 + 3.5 I, 2.6 - 3.1 I}, {5.3 - 0.5 I, 3.1 + 5.6 I}, {2.3 - 0.5 I, 3.1 + 5.6 I}};

Then the command

Covariance[Data]

gives as output the matrix

{{6.13568 + 0. I, -8.76136 - 2.69864 I}, {-8.76136 + 2.69864 I,   22.0442 + 0. I}}

when it should really give

{{6.13568 + 0. I, -8.76136 + 2.69864 I}, {-8.76136 - 2.69864 I,   22.0442 + 0. I}}

which is the output of either of the commands

n = Dimensions[Data][[1]];
1/(n - 1)*Transpose[Conjugate[Data]].(IdentityMatrix[n] - 1/n*ConstantArray[1, {12, 12}]).Data

or

1/(n - 1)*Transpose[Conjugate[(IdentityMatrix[n] - 1/n*ConstantArray[1, {12, 12}]).Data]].IdentityMatrix[n] - 1/n*ConstantArray[1, {12, 12}]).Data

What Mathematica gives is the output of

1/(n - 1)*Transpose[Data].(IdentityMatrix[n] - 1/n*ConstantArray[1, {12, 12}]).Conjugate[Data]

which is wrong. It rather seems that Mathematica is computing the Covariance matrix of the conjugate of Data.
Can someone check what is really going on?

Complex Data Problem.nb (16.9 KB)

2

Using what I think are the usual definitions of sample covariance for complex numbers, I get the same results as Mathematica:

data = {{1.2 + 0.5  I, 2.6 - 3.1  I}, {5.3 - 0.5  I, 1.1 + 5.6  I}, {3.5 + 3.5  I, 1.6 - 3.1  I}, 
            {2.5 + 3.5  I, 4.6 - 3.1  I}, {4.3 - 0.5  I, 2.1 + 5.6  I}, {3.5 + 3.5  I, 4.6 - 3.1  I}, 
            {4.3 - 0.5  I, 2.1 + 5.6  I}, {1.5 + 3.5  I, 3.6 - 3.1  I}, {2.3 - 0.5  I, 1.1 + 5.6  I}, 
            {1.5 + 3.5  I, 2.6 - 3.1  I}, {5.3 - 0.5  I, 3.1 + 5.6  I}, {2.3 - 0.5  I, 3.1 + 5.6  I}};

n = Length[data];
(* \[CapitalSigma][[1,2]] *)
(n/(n - 1)) Mean[(data[[All, 1]] - Mean[data[[All, 1]]]) Conjugate[data[[All, 2]] - Mean[data[[All, 2]]]]]
(* -8.76136 - 2.69864 I *)

(* \[CapitalSigma][2,1] *)
(n/(n - 1)) Mean[(data[[All, 2]] - Mean[data[[All, 2]]]) Conjugate[data[[All, 1]] - Mean[data[[All, 1]]]]]
(* -8.76136 + 2.69864 I *)

(* \[CapitalSigma][1,1] *)
(n/(n - 1)) Mean[(data[[All, 1]] - Mean[data[[All, 1]]]) Conjugate[data[[All, 1]] - Mean[data[[All, 1]]]]]
(* 6.13568 + 0. I *)

(* \[CapitalSigma][2,2] *)
(n/(n - 1)) Mean[(data[[All, 2]] - Mean[data[[All, 2]]]) Conjugate[data[[All, 2]] - Mean[data[[All, 2]]]]]
(* 22.0442 + 0. I *)
3

Do you have a reference for the calculation formula you use? (In other words, not in terms of Mathematica code.)