Re: covariance matrix

[Vic Norton]

On 2/20/07, at 1:40 AM -0500, John W. Eaton wrote:
> On 19-Feb-2007, Vic Norton wrote:
> | # The Octave cov function is incorrect.
> 
> OK, thanks, but I don't have time now to fix these problems myself, so
> I will have to wait for someone to submit a working patch.
> 
> Thanks,
> 
> jwe

Sorry about this, John. I went off half-cocked. There is nothing wrong with 
Octave's covariance function. I am appending a definition that I view as 
correct (accept that, in practice, I prefer to divide by m rather than m - 1 at 
the end). It agrees exactly with what Octave has.

As related to the original question, the covariance two numbers is undefined if 
you divide by m - 1 (Matlab and Octave) and 0 if you divide by m. The 
covariance matrix of row vectors of length n1 and n2 should be undefined if you 
divide by m - 1 and zeros(n1, n2) if you divide by m. However Matlab and Octave 
only allow n1 = n2 in this case; then they compute the covariance of the 
transpose of the vectors.

Regards,

Vic

#!/usr/local/bin/octave

% test1.m  - Test Octave constructs

1;

## definition of covariance function ##

function    c = cov1(x, y)
#-------------------------------------------------------------------#
#
#  Usage:   c = cov1(x, y)
#
#  Purpose: To compute the covariance of x and y
#
#-------------------------------------------------------------------#

if !(nargin == 1 || nargin == 2)
    usage("c = cov1(x)  or  c = cov1(x, y)");
endif
## note: cov1(x) = cov1(x, x).

[ m, n ] = size(x);
if m == 1, error("row dimension must be greater than 1"), endif

if nargin == 1                          # cov1(x) will be n x n
   xdev = x - ones(m, 1) * (sum(x)/m);  # deviations from the mean
   c = (xdev' * xdev)/(m - 1);          # some prefer division by m

else                                    # cov1(x, y) will be n x n1 
   [ m1, n1 ] = size(y);
   if m1 != m
      error("row dimensions do not match")
   endif
   xdev = x - ones(m, 1) * (sum(x)/m);  # deviations from the mean
   ydev = y - ones(m, 1) * (sum(y)/m);  # deviations from the mean
   c = (xdev' * ydev)/(m - 1);          # some prefer division by m
endif

endfunction



## test of covariance function ##

m = 6; n = 4;
X = rand(m, n);

printf("tests of covariance functions\n");

printf("\ntest matrix X\n");
for i = 1 : m
   printf("%10.5f", X(i, :));
   printf("\n");
endfor

printf("\ncov(X)\n");
C = cov(X);
for i = 1 : n
   printf("%10.5f", C(i, :));
   printf("\n");
endfor
printf("\ncov(X(:, 2), X(:, 3)) =%10.5f\n", cov(X(:, 2), X(:, 3)) );

printf("\ncov1(X)\n");
C = cov1(X);
for i = 1 : n
   printf("%10.5f", C(i, :));
   printf("\n");
endfor
printf("\ncov1(X(:, 2), X(:, 3)) =%10.5f\n", cov1(X(:, 2), X(:, 3)) );

n1 = 5;
Y = rand(m, n1);

printf("\ntest matrix Y\n");
for i = 1 : m
   printf("%10.5f", Y(i, :));
   printf("\n");
endfor

printf("\ncov(X, Y)\n");
C = cov(X, Y);
for i = 1 : n
   printf("%10.5f", C(i, :));
   printf("\n");
endfor

printf("\ncov1(X, Y)\n");
C = cov1(X, Y);
for i = 1 : n
   printf("%10.5f", C(i, :));
   printf("\n");
endfor