3ba3659412be25d26597e67ae1071184d9ec256e
1 // Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
2 // Copyright (C) 2013 - Scilab Enteprises - Paul Bignier: added mean squared deviation
3 //                                                        (third input argument)
4 // Copyright (C) 2000 - INRIA - Carlos Klimann
5 //
6 // Copyright (C) 2012 - 2016 - Scilab Enterprises
7 //
8 // This file is hereby licensed under the terms of the GNU GPL v2.0,
9 // pursuant to article 5.3.4 of the CeCILL v.2.1.
10 // This file was originally licensed under the terms of the CeCILL v2.1,
11 // and continues to be available under such terms.
12 // For more information, see the COPYING file which you should have received
13 // along with this program.
14 //
16 function sd = stdev(x, o, m)
17     //
18     //This function computes  the  standard deviation  of  the values  of  a
19     //vector or matrix x.
20     //
21     //For a vector or a  matrix x, s=stdev(x)  returns in  the scalar s  the
22     //standard deviation of all the entries of x.
23     //
24     //s=stdev(x,'r')  (or,  equivalently,   s=stdev(x,1))   is the   rowwise
25     //standard deviation. It returns  in each entry of the  row vector s the
26     //standard deviation of each column of x.
27     //
28     //s=stdev(x,'c')  (or,  equivalently, s=stdev(x,2))  is   the columnwise
29     //standard  deviation. It returns in each  entry of the  column vector y
30     //the standard deviation of each row of x.
31     //
32     //The input argument m represents the a priori mean. If it is present, then the sum is
33     //divided by size(x,"*"). Otherwise ("sample standard deviation"), it is divided by size(x,"*")-1.
34     //
36     [lhs, rhs] = argn(0);
38     if rhs < 1 then
39         error(msprintf(_("%s: Wrong number of input arguments: %d to %d expected.\n"),"stdev",1,3));
40     end
42     if rhs < 2 then
43         o = "*"
44         on = 0
45     else
46         select o
47         case "*" then
48             on = 0
49         case "r" then
50             on = 1
51         case "c" then
52             on = 2
53         else
54             if type(o) <> 1 | size(o, "*") <> 1 | o < 0 | floor(o) <> o then
55                 error(msprintf(_("%s: Wrong value for input argument #%d: ''%s'', ''%s'', ''%s'' or a positive integer expected.\n"),"stdev",2,"*","r","c")),
56             else
57                 on = o
58             end
59         end
60     end
62     if (length(size(x)) > 2) then
63         if rhs == 3 then
64             sd = %hm_stdev(x, o, m);
65         else
66             sd = %hm_stdev(x, o);
67         end
68         return
69     end
71     if type(x) ~= 1 | ~isreal(x) then
72         tmp = _("%s: Wrong type for input argument #%d: A real matrix expected.\n")
73         error(msprintf(tmp, "stdev", 1))
74     end
76     siz = size(x);
77     if rhs == 3 then
78         if and(typeof(m) ~= ["constant" "hypermat"]) | ~isreal(m) then
79             tmp = _("%s: Wrong type for input argument #%d: A real matrix expected.\n")
80             error(msprintf(tmp, "stdev", 3))
81         elseif on == 0 then
82             if ~isscalar(m) then
83                 tmp = _("%s: Wrong size for input argument #%d.\n")
84                 error(msprintf(tmp, "stdev", 3))
85             end
86         elseif on == 1 then
87             if or(size(m) ~= [1 siz(2:\$)]) then
88                 if isscalar(m) then
89                     m = resize_matrix(m, [1 siz(2:\$)], m);
90                 else
91                     tmp = _("%s: Wrong size for input argument #%d.\n")
92                     error(msprintf(tmp, "stdev", 3))
93                 end
94             end
95         elseif on == 2 then
96             if or(size(m) ~= [siz(1) 1 siz(3:\$)]) then
97                 if isscalar(m) then
98                     m = resize_matrix(m, [siz(1) 1 siz(3:\$)], m);
99                 else
100                     tmp = _("%s: Wrong size for input argument #%d.\n")
101                     error(msprintf(tmp, "stdev", 3))
102                 end
103             end
104         end
105     end
107     if x == [] then
108         sd = %nan;
109         return
110     end
112     if rhs == 3 & isnan(m) then
113         // This will compute the "biased variance": the denominator will be size(x,orien)
114         // but the a priori mean is not considered as provided.
115         rhs = 2
116     end
118     // Remove the mean
119     if on == 0 then
120         if rhs == 3
121             y = x - m;
122         else
123             y = x - mean(x);
124         end
125     elseif on == 2 then
126         if rhs == 3
127             y = x - m*ones(x(1,:));
128         else
129             y = x - mean(x,on)*ones(x(1,:));
130         end
131     elseif on == 1 then
132         if rhs == 3
133             y = x - ones(x(:,1))*m;
134         else
135             y = x - ones(x(:,1))*mean(x,1);
136         end
137     else // on > ndims(x) then
138         sd = 0*x;
139         return
140     end
142     mn = size(x, o);
144     if mn == 1 then
145         sd = 0*y;
146     else
147         if rhs <= 2 & exists("m", "local") == 0 then // If m is provided but rhs=2, that means we want the biased deviation
148             sd = sqrt(sum(y.^2,o)/(mn-1));
149         else
150             sd = sqrt(sum(y.^2,o)/mn);
151         end
152     end
154 endfunction