Covariance Matrix

Symbolic Computation

X1 is the first Random Variable, X2 secomd Random Variable and so on the sequence X1 X2 X3 ... a sequence of Random Variables of any length. 

X1j for a fixed column j is the value of Xj happening with the value X1j at first trial, X2j the value of Xj happening at 2nd trial and so on. 

The Random Variables X1, X2, X3 ... form a 2 dimensional matrix:

1. Each column is the set of repeated values of the Random Variable given some probability ditribution 
2. Each row is a sample of an event which is happening of X1 followed by happening of X2 and so on.

(* Treat the Random Variable sequence as a matrix of numbers *)
X = {{X11, X12, X13}, {X21, X22, X23}, {X31, X32, X33}};
TableForm[X,TableHeadings->{{"Trial 1","Trial 2", "Trial 3"},{"X1","X2","X3"}}]
X//MatrixForm
Mean[X]//MatrixForm

(* Find the covariance of the first column and second column, treating the Random Variables as number sequences  *)
covx12=Covariance[{X11, X21, X31}, {X12, X22, X32}];
Simplify[covx12,Element[{X12, X22, X32},Reals]]

Or

(* X[[All, 1]] means column 1 of X *)
covx12 =Covariance[X[[All, 1]], X[[All, 2]]];
Simplify[covx12,Element[{X12, X22, X32},Reals]]

(* Treat the Random Variables as a matrix, and calculate variances of all possible columns *)
covx=Covariance[X]; 
(* Show the covariance of first and second column *)
Simplify[covx[[1]][[2]],Element[{X12, X22, X32},Reals]]

Let’s compute the covariance of two matrices against each other:

Y = {{Y11, Y12, Y13}, {Y21, Y22, Y23}, {Y31, Y32, Y33}};
Y//MatrixForm
covxy = Covariance [X, Y];

Covariance matrix is index by ij which hold the covariance of ith column of X against jth column of Y:

i = 2;
j = 3;
covxyij=Covariance [X[[All, i]], Y[[All, j]]];
Expand[Simplify[covxyij,Element[{Y22, Y13, Y33, Y23},Reals]]]
Expand[Simplify[covxy[[i,j]],Element[{Y22, Y13, Y33, Y23},Reals]]]

There is another way of calculating the entries of this matrix:

Simplify[Mean[X[[All,i]]*Y[[All,j]]]-Mean[X[[All,i]]*Mean[Y[[All,j]]]]- (covxy[[i,j]]*((Length[X[[All,i]]]-1)/Length[X[[All,i]]])),Element[{Y22, Y13, Y33, Y23},Reals]]

Remark: In the book the covariance divides by 1/n in Mathematica it divides by 1/(n-1) so adjustment was added.

Numeric Computation

n = 9;
(* X is (X1, X2, X3, ..., Xn)' p.49 *)
X = Transpose[Table[RandomVariate[NormalDistribution[1,3],n], {i, 1, n}]];
TableForm[X,TableHeadings->{{"Trial 1","Trial 2", "Trial 3","Trial 4","Trial 5", "Trial 6", "Trial 7","Trial 8", "Trial 9"},{"X1","X2","X3","X4","X5","X6","X7","X8","X9"}}]
Mean[X]
covx=Covariance[X];
covx//MatrixForm

2D Visualization 

Increase n to a number larger than 30, the diagonal is distinguished for higher values:

n = 10^2;
X = Transpose[Table[RandomVariate[NormalDistribution[1,3],n], {i, 1, n}]];
covx=Covariance[X];
covx//MatrixForm;
MatrixPlot[covx]
ArrayPlot[covx]

3D Visualization

n = 10^3;
m = 1000; (* X need not be a square matrix as above, since it will be computationally too expensive, so we can choose much smaller numbers for the horizontal axis i.e. 1, 2, ..., m << n *)
σ=4;
μ = 1;
X = Transpose[Table[RandomVariate[NormalDistribution[μ,σ],m], {i, 1, n}]];
(* Notice the mean is around 1 as it was passed to the NormalDistribution in the previous line *)
mean=Mean[X];
ListPlot[mean]
covx=Covariance[X];
(* covx convariance matrix is 1000x1000 matrix of reals, no way to display, so let's display only first 30 rows and columns *)
(* Notice the diagonal of the matrix is Covariance(X, X) which is Var(X) which is σ *)
ListPlot3D[covx[[1;;30, 1;;30]],Mesh->None,InterpolationOrder->0,ColorFunction->"SouthwestColors", ImageSize->500, PlotRange->All]

Reduce the n, m ( < 30 ) the number of samples and you can observe that the diagonal is no longer distinguished from the rest of the matrix:

n = 20;
m = 20; 
σ=4;
μ = 1;
X = Transpose[Table[RandomVariate[NormalDistribution[μ,σ],m], {i, 1, n}]];
mean=Mean[X];
ListPlot[mean]
covx=Covariance[X];
ListPlot3D[covx[[1;;20, 1;;20]],Mesh->None,InterpolationOrder->0,ColorFunction->"SouthwestColors", ImageSize->500, PlotRange->All]

As n, m larger the diagonal is more pronounced:

n = 50;
m = 50; 
σ=4;
μ = 1;
X = Table[RandomVariate[NormalDistribution[μ,σ],m], {i, 1, n}];
mean=Mean[X];
ListPlot[mean]
covx=Covariance[X];
ListPlot3D[covx[[1;;30, 1;;30]],Mesh->None,InterpolationOrder->0,ColorFunction->"SouthwestColors", ImageSize->500, PlotRange->All]