/*  -*- MaTX -*-
 *
 *  NAME
 *      cov()     - Covariance
 *      cov_col() - Covariance matrix for each column
 *      cov_row() - Covariance matrix for each row
 *
 *  SYNOPSIS
 *      cv = cov(x)
 *          Real cv;
 *          Matrix x;
 *
 *      cv = cov(x, y)
 *          Real cv;
 *          Matrix x;
 *          Matrix y;
 *
 *      cv = cov_col(x)
 *          Matrix cv;
 *          Matrix x;
 *
 *      cv = cov_col(x, y)
 *          Matrix cv;
 *          Matrix x;
 *          Matrix y;
 *
 *      cv = cov_row(x)
 *          Matrix cv;
 *          Matrix x;
 *
 *      cv = cov_row(x, y)
 *          Matrix cv;
 *          Matrix x;
 *          Matrix y;
 *
 *  DESCRIPTION
 *      cov(x) returns the variance for the elements of x.
 *      cov(x,y) is equal to cov([x y]).  
 *
 *      If x is a matrix whose each row is an observation, and each 
 *      column a variable, cov_col(x) is the covariance matrix. 
 *
 *      diag_vec(cov_col(x)) is a vector of variances for each column,
 *      and sqrt(diag_vec(cov_col(x))) is the standard deviations.
 *
 *      If x is a matrix whose each column is an observation, and each 
 *      row a variable, cov_row(x) is the covariance matrix. 
 *
 *      diag_vec(cov_row(x)) is a vector of variances for each row,
 *      and sqrt(diag_vec(cov_row(x))) is the standard deviations.
 *
 *  SEE ALSO
 *      corrcoef and std
 */

Func Real cov(x, y, ...)
	Matrix x, y;
{
	Matrix z;

	error(nargchk(1, 2, nargs, "cov"));

	if (nargs == 1) {
		return [cov_col(reshape(x, 1, Rows(x)*Cols(x)))](1,1);
	} else {
		z = [makecolv(x), makecolv(y)];
		return [cov_col(reshape(z, 1, Rows(z)*Cols(z)))](1,1);
	}
}

Func Matrix cov_col(x_, y_, ...)
	Matrix x_, y_;
{
	Integer m, n;
	Matrix X, cv;

	error(nargchk(1, 2, nargs, "cov_col"));

	{m, n} = size(x_);

	if (nargs > 1) {
		if (m != Rows(y_) || n != Cols(y_)) {
			error("cov_col(): two vectors must be the same size.\n");
		}
		if (m != 1 && n != 1) {
			error("cov_col(): arguments must be vectors.\n");
		}
		X = [makecolv(x_), makecolv(y_)];
	} else if (min(Cols(x_), Rows(x_)) == 1) {
		X = makecolv(x_);
	} else {
		X = x_;
	}

	m = Rows(X);
	X = X - ONE(m,1)*(sum_col(X)/m);

	if (m == 1) {
		cv = Z(1);
	} else {
		cv = X#*X / (m-1);
	}
	return cv;
}

Func Matrix cov_row(x, y, ...)
	Matrix x, y;
{
	error(nargchk(1, 2, nargs, "cov_row"));
	
	if (nargs == 1) {
		return trans(cov_col(trans(x)));
	} else {
		return trans(cov_col(x, y));
	}
}