* Bug #13050 fixed - The result of mvvacov function was not symmetric.
[scilab.git] / scilab / modules / statistics / macros / cov.sci
1 // Copyright (C) 2012-2013 - Michael Baudin
2 // Copyright (C) 2009-2010 - DIGITEO - Michael Baudin
3 // Copyright (C) 1993 - 1995 - Anders Holtsberg
4 //
5 // This file must be used under the terms of the CeCILL.
6 // This source file is licensed as described in the file COPYING, which
7 // you should have received as part of this distribution.  The terms
8 // are also available at
9 // http://www.cecill.info/licences/Licence_CeCILL_V2.1-en.txt
10
11 function C = cov(varargin)
12     // Covariance matrix
13     //
14     // Calling Sequence
15     //   C = cov(x)
16     //   C = cov(x, 0)
17     //   C = cov(x, 1)
18     //   C = cov(x, y)
19     //   C = cov(x, y, 0)
20     //   C = cov(x, y, 1)
21     //
22     // Parameters
23     // x: a nobs-by-1 or nobs-by-nvar matrix of doubles
24     // y: a nobs-by-1 or nobs-by-nvar matrix of doubles
25     // C: a square matrix of doubles, the empirical covariance
26     //
27     // Description
28     // If x is a nobs-by-1 matrix,
29     // then cov(x) returns the variance of x,
30     // normalized by nobs-1.
31     //
32     // If x is a nobs-by-nvar matrix,
33     // then cov(x) returns the nvar-by-nvar covariance matrix of the
34     // columns of x, normalized by nobs-1.
35     // Here, each column of x is a variable and
36     // each row of x is an observation.
37     //
38     // If x and y are two nobs-by-1 matrices,
39     // then cov(x, y) returns the 2-by-2 covariance matrix of x and
40     // y, normalized by nobs-1, where nobs is the number of observations.
41     //
42     // cov(x, 0) is the same as cov(x) and
43     // cov(x, y, 0) is the same as cov(x, y).
44     // In this case, if the population is from a normal distribution,
45     // then C is the best unbiased estimate of the covariance matrix.
46     //
47     // cov(x, 1) and cov(x, y, 1) normalize by nobs.
48     // In this case, C is the second moment matrix of the
49     // observations about their mean.
50     //
51     // The covariance of X and Y is defined by:
52     //
53     // Cov(X, Y) = E( (X-E(X)) (Y-E(Y))^T )
54     //
55     // where E is the expectation.
56     //
57     // This function is compatible with Matlab.
58     //
59     // Examples
60     // x = [1; 2];
61     // y = [3; 4];
62     // C = cov(x, y)
63     // expected = [0.5, 0.5; 0.5, 0.5]
64     // C = cov([x, y])
65     //
66     // x = [230; 181; 165; 150; 97; 192; 181; 189; 172; 170];
67     // y = [125; 99; 97; 115; 120; 100; 80; 90; 95; 125];
68     // expected = [
69     // 1152.4556, -88.911111
70     // -88.911111, 244.26667
71     // ]
72     // C = cov(x, y)
73     // C = cov([x, y])
74     //
75     // // Source [3]
76     // A = [
77     // 4.0 2.0 0.60
78     // 4.2 2.1 0.59
79     // 3.9 2.0 0.58
80     // 4.3 2.1 0.62
81     // 4.1 2.2 0.63
82     // ];
83     // S = [
84     // 0.025 0.0075 0.00175
85     // 0.0075 0.007 0.00135
86     // 0.00175 0.00135 0.00043
87     // ];
88     // C = cov(A)
89     //
90     // Bibliography
91     // [1] http://en.wikipedia.org/wiki/Covariance_matrix
92     // [2] "Introduction to probability and statistics for engineers and scientists.", Sheldon Ross
93     // [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
94
95
96     [lhs, rhs]=argn()
97     //
98     if (rhs == 1) then
99         x = varargin(1)
100         //
101         // Check type
102         if (typeof(x) <> "constant")
103             error(msprintf(gettext("%s: Wrong type for input argument #%d: a real matrix expected.\n"),"cov", 1));
104         end
105         nobs = size(x, "r")
106         r = 1/(nobs-1)
107         A = x
108     elseif (rhs == 2) then
109         //
110         x = varargin(1)
111         y = varargin(2)
112         //
113         // Check type
114         if (typeof(x) <> "constant")
115             error(msprintf(gettext("%s: Wrong type for input argument #%d: a real matrix expected.\n"),"cov", 1));
116         end
117         if (typeof(y) <> "constant")
118             error(msprintf(gettext("%s: Wrong type for input argument #%d: an integer or a real matrix expected.\n"),"cov", 2));
119         end
120         //
121         // Check size
122         nobs = size(x, "r")
123         if (size(y, "*") == 1) then
124             if (y <> 0 & y <> 1)
125                 error(msprintf(gettext("%s: Wrong value for input argument #%d: %d or %d expected.\n"),"cov", 2, 0, 1));
126             end
127             if (y == 1) then
128                 r = 1/nobs
129                 A = x
130             elseif (y == 0) then
131                 r = 1/(nobs-1)
132                 A = x
133             end
134         else
135             if (size(x) <> [nobs 1]) then
136                 error(msprintf(gettext("%s: Wrong size for input argument #%d: %dx%d expected.\n"),"cov", 1, nobs, 1));
137             end
138             if (size(y) <> [nobs 1]) then
139                 error(msprintf(gettext("%s: Wrong size for input argument #%d: %dx%d expected.\n"),"cov", 2, nobs, 1));
140             end
141             r = 1/(nobs-1)
142             A = [x, y]
143         end
144     elseif (rhs == 3) then
145         //
146         x = varargin(1)
147         y = varargin(2)
148         nrmlztn = varargin(3)
149         //
150         // Check type
151         if (typeof(x) <> "constant")
152             error(msprintf(gettext("%s: Wrong type for input argument #%d: a real matrix expected.\n"),"cov", 1));
153         end
154         if (typeof(y) <> "constant")
155             error(msprintf(gettext("%s: Wrong type for input argument #%d: a real matrix expected.\n"),"cov", 2));
156         end
157         if (typeof(nrmlztn) <> "constant")
158             error(msprintf(gettext("%s: Wrong type for input argument #%d: an integer expected.\n"),"cov", 3));
159         end
160         //
161         // Check size
162         nobs = size(x, "r")
163         if (size(x) <> [nobs 1]) then
164             error(msprintf(gettext("%s: Wrong size for input argument #%d: %dx%d expected.\n"),"cov", 1, nobs, 1));
165         end
166         if (size(y) <> [nobs 1]) then
167             error(msprintf(gettext("%s: Wrong size for input argument #%d: %dx%d expected.\n"),"cov", 2, nobs, 1));
168         end
169         if (size(nrmlztn, "*") <> 1) then
170             error(msprintf(gettext("%s: Wrong type for input argument #%d: an integer expected.\n"),"cov", 3));
171         end
172         //
173         // Check content
174         if (nrmlztn <> 0 & nrmlztn <> 1)
175             error(msprintf(gettext("%s: Wrong value for input argument #%d: %d or %d expected.\n"),"cov", 3, 0, 1));
176         end
177         A = [x, y]
178         if (nrmlztn == 1) then
179             r = 1/nobs
180         else
181             r = 1/(nobs-1)
182         end
183     else
184         error(msprintf(gettext("%s: Wrong number of input argument(s): %d, %d or %d expected.\n"),"cov", 1, 2, 3));
185     end
186     //
187     // Compute with A in the general case
188     nvar = size(A, "c")
189     nobs = size(A, "r")
190     for i = 1:nvar
191         A(:,i) = A(:,i) - mean(A(:,i))
192     end
193     C = zeros(nvar, nvar)
194     for i = 1:nvar
195         C(i,i) = A(:,i)'*A(:,i)*r
196         for j = i+1:nvar
197             C(i,j) = A(:,i)'*A(:,j)*r
198             C(j,i) = C(i,j)
199         end
200     end
201 endfunction