[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]
Re: covariance matrix
From: |
Vic Norton |
Subject: |
Re: covariance matrix |
Date: |
Tue, 20 Feb 2007 10:27:41 -0500 |
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