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         msg = _("%s: Wrong number of input arguments: %d to %d expected.\n")
40         error(msprintf(msg, "stdev", 1, 3))
41     end
43     if rhs < 2 then
44         o = "*"
45         on = 0
46     else
47         select o
48         case "*" then
49             on = 0
50         case "r" then
51             on = 1
52         case "c" then
53             on = 2
54         else
55             if type(o) <> 1 || size(o, "*") <> 1 || floor(o) <> o | o < 1 | o > ndims(x) then
56                 msg = _("%s: Argument #%d: Must be in the set {%s}.\n")
57                 sset =  """*"" ""r"" ""c"" 1 2"
58                 if ndims(x)>2
59                     sset = sset + strcat(msprintf(" %d\n",(3:ndims(x))'))
60                 end
61                 error(msprintf(msg, "stdev", 2, sset))
62             else
63                 on = o
64             end
65         end
66     end
68     if (length(size(x)) > 2) then
69         if rhs == 3 then
70             sd = %hm_stdev(x, on, m);
71         else
72             sd = %hm_stdev(x, on);
73         end
74         return
75     end
77     if type(x) ~= 1 | ~isreal(x) then
78         tmp = _("%s: Wrong type for input argument #%d: A real matrix expected.\n")
79         error(msprintf(tmp, "stdev", 1))
80     end
82     siz = size(x);
83     if rhs == 3 then
84         if typeof(m) ~= "constant" | ~isreal(m) then
85             tmp = _("%s: Wrong type for input argument #%d: A real matrix expected.\n")
86             error(msprintf(tmp, "stdev", 3))
87         elseif on == 0 then
88             if ~isscalar(m) then
89                 tmp = _("%s: Wrong size for input argument #%d.\n")
90                 error(msprintf(tmp, "stdev", 3))
91             end
92         elseif on == 1 then
93             if or(size(m) ~= [1 siz(2:\$)]) then
94                 if isscalar(m) then
95                     m = resize_matrix(m, [1 siz(2:\$)], m);
96                 else
97                     tmp = _("%s: Wrong size for input argument #%d.\n")
98                     error(msprintf(tmp, "stdev", 3))
99                 end
100             end
101         elseif on == 2 then
102             if or(size(m) ~= [siz(1) 1 siz(3:\$)]) then
103                 if isscalar(m) then
104                     m = resize_matrix(m, [siz(1) 1 siz(3:\$)], m);
105                 else
106                     tmp = _("%s: Wrong size for input argument #%d.\n")
107                     error(msprintf(tmp, "stdev", 3))
108                 end
109             end
110         end
111     end
113     if x == [] then
114         sd = %nan;
115         return
116     end
118     if rhs == 3 & isnan(m) then
119         // This will compute the "biased variance": the denominator will be size(x,orien)
120         // but the a priori mean is not considered as provided.
121         rhs = 2
122     end
124     // Remove the mean
125     if on == 0 then
126         if rhs == 3
127             y = x - m;
128         else
129             y = x - mean(x);
130         end
131     elseif on == 2 then
132         if rhs == 3
133             y = x - m*ones(x(1,:));
134         else
135             y = x - mean(x,on)*ones(x(1,:));
136         end
137     elseif on == 1 then
138         if rhs == 3
139             y = x - ones(x(:,1))*m;
140         else
141             y = x - ones(x(:,1))*mean(x,1);
142         end
143     else // on > ndims(x) then
144         sd = 0*x;
145         return
146     end
148     mn = size(x, o);
150     if mn == 1 then
151         sd = 0*y;
152     else
153         if rhs <= 2 & exists("m", "local") == 0 then // If m is provided but rhs=2, that means we want the biased deviation
154             sd = sqrt(sum(y.^2,o)/(mn-1));
155         else
156             sd = sqrt(sum(y.^2,o)/mn);
157         end
158     end
160 endfunction