// Copyright (C) 2012-2013 - Michael Baudin
// Copyright (C) 2009-2010 - DIGITEO - Michael Baudin
// Copyright (C) 1993 - 1995 - Anders Holtsberg
//
// This file must be used under the terms of the CeCILL.
// This source file is licensed as described in the file COPYING, which
// you should have received as part of this distribution. The terms
// are also available at
// http://www.cecill.info/licences/Licence_CeCILL_V2.1-en.txt
function C = cov(varargin)
// Covariance matrix
//
// Calling Sequence
// C = cov(x)
// C = cov(x, 0)
// C = cov(x, 1)
// C = cov(x, y)
// C = cov(x, y, 0)
// C = cov(x, y, 1)
//
// Parameters
// x: a nobs-by-1 or nobs-by-nvar matrix of doubles
// y: a nobs-by-1 or nobs-by-nvar matrix of doubles
// C: a square matrix of doubles, the empirical covariance
//
// Description
// If x is a nobs-by-1 matrix,
// then cov(x) returns the variance of x,
// normalized by nobs-1.
//
// If x is a nobs-by-nvar matrix,
// then cov(x) returns the nvar-by-nvar covariance matrix of the
// columns of x, normalized by nobs-1.
// Here, each column of x is a variable and
// each row of x is an observation.
//
// If x and y are two nobs-by-1 matrices,
// then cov(x, y) returns the 2-by-2 covariance matrix of x and
// y, normalized by nobs-1, where nobs is the number of observations.
//
// cov(x, 0) is the same as cov(x) and
// cov(x, y, 0) is the same as cov(x, y).
// In this case, if the population is from a normal distribution,
// then C is the best unbiased estimate of the covariance matrix.
//
// cov(x, 1) and cov(x, y, 1) normalize by nobs.
// In this case, C is the second moment matrix of the
// observations about their mean.
//
// The covariance of X and Y is defined by:
//
// Cov(X, Y) = E( (X-E(X)) (Y-E(Y))^T )
//
// where E is the expectation.
//
// This function is compatible with Matlab.
//
// Examples
// x = [1; 2];
// y = [3; 4];
// C = cov(x, y)
// expected = [0.5, 0.5; 0.5, 0.5]
// C = cov([x, y])
//
// x = [230; 181; 165; 150; 97; 192; 181; 189; 172; 170];
// y = [125; 99; 97; 115; 120; 100; 80; 90; 95; 125];
// expected = [
// 1152.4556, -88.911111
// -88.911111, 244.26667
// ]
// C = cov(x, y)
// C = cov([x, y])
//
// // Source [3]
// A = [
// 4.0 2.0 0.60
// 4.2 2.1 0.59
// 3.9 2.0 0.58
// 4.3 2.1 0.62
// 4.1 2.2 0.63
// ];
// S = [
// 0.025 0.0075 0.00175
// 0.0075 0.007 0.00135
// 0.00175 0.00135 0.00043
// ];
// C = cov(A)
//
// Bibliography
// [1] http://en.wikipedia.org/wiki/Covariance_matrix
// [2] "Introduction to probability and statistics for engineers and scientists.", Sheldon Ross
// [3] NIST/SEMATECH e-Handbook of Statistical Methods, 6.5.4.1. Mean Vector and Covariance Matrix, http://www.itl.nist.gov/div898/handbook/pmc/section5/pmc541.htm
[lhs, rhs]=argn()
//
if (rhs == 1) then
x = varargin(1)
//
// Check type
if (typeof(x) <> "constant")
error(msprintf(gettext("%s: Wrong type for input argument #%d: a real matrix expected.\n"),"cov", 1));
end
nobs = size(x, "r")
r = 1/(nobs-1)
A = x
elseif (rhs == 2) then
//
x = varargin(1)
y = varargin(2)
//
// Check type
if (typeof(x) <> "constant")
error(msprintf(gettext("%s: Wrong type for input argument #%d: a real matrix expected.\n"),"cov", 1));
end
if (typeof(y) <> "constant")
error(msprintf(gettext("%s: Wrong type for input argument #%d: an integer or a real matrix expected.\n"),"cov", 2));
end
//
// Check size
nobs = size(x, "r")
if (size(y, "*") == 1) then
if (y <> 0 & y <> 1)
error(msprintf(gettext("%s: Wrong value for input argument #%d: %d or %d expected.\n"),"cov", 2, 0, 1));
end
if (y == 1) then
r = 1/nobs
A = x
elseif (y == 0) then
r = 1/(nobs-1)
A = x
end
else
if (size(x) <> [nobs 1]) then
error(msprintf(gettext("%s: Wrong size for input argument #%d: %dx%d expected.\n"),"cov", 1, nobs, 1));
end
if (size(y) <> [nobs 1]) then
error(msprintf(gettext("%s: Wrong size for input argument #%d: %dx%d expected.\n"),"cov", 2, nobs, 1));
end
r = 1/(nobs-1)
A = [x, y]
end
elseif (rhs == 3) then
//
x = varargin(1)
y = varargin(2)
nrmlztn = varargin(3)
//
// Check type
if (typeof(x) <> "constant")
error(msprintf(gettext("%s: Wrong type for input argument #%d: a real matrix expected.\n"),"cov", 1));
end
if (typeof(y) <> "constant")
error(msprintf(gettext("%s: Wrong type for input argument #%d: a real matrix expected.\n"),"cov", 2));
end
if (typeof(nrmlztn) <> "constant")
error(msprintf(gettext("%s: Wrong type for input argument #%d: an integer expected.\n"),"cov", 3));
end
//
// Check size
nobs = size(x, "r")
if (size(x) <> [nobs 1]) then
error(msprintf(gettext("%s: Wrong size for input argument #%d: %dx%d expected.\n"),"cov", 1, nobs, 1));
end
if (size(y) <> [nobs 1]) then
error(msprintf(gettext("%s: Wrong size for input argument #%d: %dx%d expected.\n"),"cov", 2, nobs, 1));
end
if (size(nrmlztn, "*") <> 1) then
error(msprintf(gettext("%s: Wrong type for input argument #%d: an integer expected.\n"),"cov", 3));
end
//
// Check content
if (nrmlztn <> 0 & nrmlztn <> 1)
error(msprintf(gettext("%s: Wrong value for input argument #%d: %d or %d expected.\n"),"cov", 3, 0, 1));
end
A = [x, y]
if (nrmlztn == 1) then
r = 1/nobs
else
r = 1/(nobs-1)
end
else
error(msprintf(gettext("%s: Wrong number of input argument(s): %d, %d or %d expected.\n"),"cov", 1, 2, 3));
end
//
// Compute with A in the general case
nvar = size(A, "c")
nobs = size(A, "r")
for i = 1:nvar
A(:,i) = A(:,i) - mean(A(:,i))
end
C = zeros(nvar, nvar)
for i = 1:nvar
C(i,i) = A(:,i)'*A(:,i)*r
for j = i+1:nvar
C(i,j) = A(:,i)'*A(:,j)*r
C(j,i) = C(i,j)
end
end
endfunction