scilab/modules/statistics/macros/variance.sci

   1 // Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
   2 // Copyright (C) 2000 - INRIA - Carlos Klimann
   3 // Copyright (C) 2013 - Samuel GOUGEON
   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
  12 function [s, m] = variance(x, orien, m)
  13     //
  14     //This function computes  the  variance  of  the values of  a  vector or
  15     //matrix x.
  16     //
  17     //For a vector or   a matrix x,  s=variance(x)  returns in  the scalar s  the
  18     //variance of all the entries of x.
  19     //
  20     //s=variance(x,'r')    (or,  equivalently,    s=variance(x,1)) is     the  rowwise
  21     //variance. It returns in each entry of the row vector s the variance of
  22     //each column of x.
  23     //
  24     //s=variance(x,'c') (or, equivalently, s=variance(x,2)) is the columnwise standard
  25     //deviation.   It  returns in  each  entry  of the   column vector y the
  26     //variance of each row of x.
  27     //
  28     //The input argument m represents the a priori mean. If it is present, then the sum is
  29     //divided by n. Otherwise ("sample variance"), it is divided by n-1.
  30     //
  31
  32     // Checking and normalizing input arguments:
  33     // ----------------------------------------
  34     [lhs,rhs] = argn(0)
  35     if rhs==0 then
  36         tmp = gettext("%s: Wrong number of input arguments: %d to %d expected.\n")
  37         error(msprintf(tmp, "variance", 1, 2))
  38     end
  39
  40     if x==[] then
  41         s = %nan
  42         return
  43     end
  44
  45     if ~isdef("orien","local") then
  46         orien = "*"
  47     end
  48
  49     if rhs==3 then
  50         if typeof(m)~="constant" then
  51             tmp = gettext("%s: Wrong value of m : a priori mean expected.\n")
  52             error(msprintf(tmp, "variance"))
  53         elseif orien=="*" then
  54             if ~isscalar(m) then
  55                 tmp = gettext("%s: Wrong value of m : a priori mean expected.\n")
  56                 error(msprintf(tmp, "variance"))
  57             end
  58         elseif orien=="r" | orien==1 then
  59             if size(m)~=[1 size(x,"c")] & ~isscalar(m) then
  60                 tmp = gettext("%s: Wrong value of m : a priori mean expected.\n")
  61                 error(msprintf(tmp, "variance"))
  62             end
  63         elseif orien=="c" | orien==2 then
  64             if size(m)~=[size(x,"r") 1] & ~isscalar(m) then
  65                 tmp = gettext("%s: Wrong value of m : a priori mean expected.\n")
  66                 error(msprintf(tmp, "variance"))
  67             end
  68         end
  69     end
  70
  71     transposed = %f // to refer and process as in "r", we priorly transpose any "c" request
  72     if orien=="r" | orien==1 | orien=="c" | orien==2 | orien=="*"
  73         if orien=="c" | orien==2 then
  74             x = x.'
  75             transposed = %t
  76             orien = "r"
  77         end
  78     else
  79         tmp = gettext("%s: Wrong value for input argument #%d: ''%s'', ''%s'', %d or %d expected.\n")
  80         error(msprintf(tmp, "variance", 2, "c", "r", 1, 2))
  81     end
  82
  83     // Calculations
  84     // ------------
  85
  86     d = size(x, orien) - 1 + exists("m","local") // Denominator. If m is given, then the a priori mean is known and we divide by size(n,orien)
  87
  88     if rhs == 3 & isnan(m) then
  89         // This will compute the "biased variance": the denominator is size(x,orien) but the a priori mean is not considered as provided.
  90         rhs = 2
  91     end
  92     if orien=="*" then
  93         if rhs < 3 then
  94             m = mean(x)
  95         end
  96     else
  97         if rhs < 3 then
  98             m = mean(x, orien).*.ones(size(x,1),1)
  99         else
 100             if isscalar(m) then
 101                 if or(m==[0 1]) then
 102                     tmp = _("%s: The significance of input argument #%d has been modified. Please refer to the variance help page.\n")
 103                     warning(msprintf(tmp, "variance", 3))
 104                 end
 105                 // If m is a scalar, extend it to the size of x.
 106                 // If lhs==1, we do not need to perform this operation, because in the following 'x - m', m can be a scalar
 107                 m = m*ones(x)
 108             else
 109                 if transposed then
 110                     m = m.';
 111                 end
 112                 m = m.*.ones(size(x,1),1)
 113             end
 114         end
 115     end
 116
 117     s = sum(abs(x - m).^2, orien) / d
 118
 119     m = m(1, :);
 120     if transposed then
 121         s = s.'
 122         m = m.'
 123     end
 124
 125 endfunction