Each `xi`

should be a vector (random variable) with it's own variance and mean.

Covariance matrix is symmetric, so you just need to compute one half of it (and copy the rest) and has variance of xi at main diagonal.

```
S = ...// your symmetric matrix n*n
for(int i=0; i<n;i++)
S(i,i) = var(xi);
for(j = i+1; j<n; j++)
S(i,j) = cov(xi, xj);
S(j,i) = S(i,j);
end
end
```

where variance (var) of xi:

```
v = 0;
for(int i = 0; i<xi.Count; i++)
v += (xi(i) - mean(xi))^2;
end
v = v / xi.Count;
```

and covariance (cov)

```
cov(xi, xj) = r(xi,xj) * sqrt(var(xi)) * sqrt(var(xj))
```

where `r(xi, xj)`

is Pearson product-moment correlation coefficient

**EDIT**

or, since cov(X, Y) = E(X*Y) - E(X)*E(Y)

```
cov(xi, xj) = mean(xi.*xj) - mean(xi)*mean(xj);
```

where `.*`

is Matlab-like element-wise multiplication.

So if `x = [x1, x2]`

, `y = [y1, y2]`

then `z = x.*y = [x1*y1, x2*y2]`

;