function [covXY] = covariance(x, y, s)
%covariance Computes a variance-covariance matrix.
%
% Computes the variance-covariance matrix.
%
% USE:
% [covVec] = covariance(x, y, s),
%
% where x = n by m matrix,
% y = n by 1 vector,
% covVec = m by 1 vector or covariance between each column of x and y.
% s = the denominator for the calculation of variance and covariance.
%     set = -1 to use n-1.  Default is n.
%
% If only one input argument was given,
% it returns the variance of x.  Variance and covariances are
% calculated by using the sample size (n).  To compute with n-1, use
% the option 1.
% 
% See also corrcoef, cov, and cor.

% 23 January 2002
% Tomo Eguchi


if nargin == 1,
   covXY = var(x, 1);        % if only one input, return variance of x.
elseif nargin > 3,
   error('Too many input arguments.');
elseif nargin == 2,
   if (length(y)==1),
      s = y;
      if s == 0,
         covXY = var(x, 1);
      elseif s == 1,
         covXY = var(x);
      end
      
   else
      s = 0;
   end
   
else
   if (s < 0 | s > 1),
      error('The third input should be either 0 or 1.');
   elseif (s > 0 & s < 1),
      error('The third input should be an integer (either 0 or 1).');
   end
end

[nrX, ncX] = size(x);
[nrY, ncY] = size(y);
den = nrX - s;

if (nrX ~= nrY),
   error('The number of observations in x is not equal to that of y.'),
end

if ncY == 1,
   X2 = x - repmat(mean(x), nrX, 1);     % subtract column means
   y2 = y - mean(y);
   Y2 = repmat(y2, 1, ncX);
   covXY = (X2'*Y2)./den;
else
   error('Too many columns in the second input.');
end