# HG changeset patch
# User Stefano Guidoni <ilguido@users.sf.net>
# Date 1630515079 -7200
#      Wed Sep 01 18:51:19 2021 +0200
# Node ID 7c284aaf814cd69c2dc36fa56b61bc0da5609c33
# Parent  8f1a480c1c6cbc9e05afb9fc708dd6a7ab5c4a13
std: weighted standard deviation

diff -r 8f1a480c1c6c -r 7c284aaf814c scripts/statistics/std.m
--- a/scripts/statistics/std.m	Wed Sep 01 18:17:30 2021 +0200
+++ b/scripts/statistics/std.m	Wed Sep 01 18:51:19 2021 +0200
@@ -25,8 +25,8 @@
 
 ## -*- texinfo -*-
 ## @deftypefn  {} {} std (@var{x})
-## @deftypefnx {} {} std (@var{x}, @var{opt})
-## @deftypefnx {} {} std (@var{x}, @var{opt}, @var{dim})
+## @deftypefnx {} {} std (@var{x}, @var{w})
+## @deftypefnx {} {} std (@var{x}, @var{w}, @var{dim})
 ## Compute the standard deviation of the elements of the vector @var{x}.
 ##
 ## The standard deviation is defined as
@@ -51,7 +51,7 @@
 ## If @var{x} is a matrix, compute the standard deviation for each column and
 ## return them in a row vector.
 ##
-## The argument @var{opt} determines the type of normalization to use.
+## The argument @var{w} determines the weighting scheme to use.
 ## Valid values are
 ##
 ## @table @asis
@@ -62,13 +62,18 @@
 ## @item 1:
 ##   normalize with @math{N}, this provides the square root of the second
 ## moment around the mean
+##
+## @item a vector:
+##   compute the weighted standard deviation with nonnegative scalar weights,
+## the length of @var{w} must be equal to the size of @var{x} along dimension
+## @var{dim}
 ## @end table
 ##
 ## If the optional argument @var{dim} is given, operate along this dimension.
 ## @seealso{var, bounds, mad, range, iqr, mean, median}
 ## @end deftypefn
 
-function retval = std (x, opt = 0, dim)
+function retval = std (x, w = 0, dim)
 
   if (nargin < 1)
     print_usage ();
@@ -78,10 +83,13 @@
     error ("std: X must be a numeric vector or matrix");
   endif
 
-  if (isempty (opt))
-    opt = 0;
-  elseif (! isscalar (opt) || (opt != 0 && opt != 1))
-    error ("std: normalization OPT must be 0 or 1");
+  if (isempty (w))
+    w = 0;
+  elseif (! isvector (w) ||
+          ! isnumeric (w) ||
+          (isvector (w) && any (w < 0)) ||
+          (isscalar (w) && (w != 0 && w != 1)))
+    error ("std: W must be 0 or 1 or a vector of nonnegative scalars");
   endif
 
   nd = ndims (x);
@@ -103,7 +111,26 @@
       retval = zeros (sz);
     endif
   else
-    retval = sqrt (sumsq (center (x, dim), dim) / (n - 1 + opt));
+    if (isscalar (w))
+      retval = sqrt (sumsq (center (x, dim), dim) / (n - 1 + w));
+    else
+      ## Weighted standard deviation 
+      if (length (w) != n)
+        error (["std: the length of W must be equal to the size of X "... 
+                "along which the standard deviation is calculated"]);
+      else
+        if ((dim == 1 && rows (w) == 1) ||
+            (dim == 2 && columns (w) == 1))
+          w = transpose (w);
+        elseif (dim > 2)
+          newdims = [(ones (1, (dim - 1))), (length (w))];
+          w = reshape (w, newdims);
+        endif
+        den = sum (w);
+        mu = sum (w .* x) ./ sum (w);
+        retval = sqrt (sum (w .* ((x .- mu) .^ 2)) / den);
+      endif
+    endif  
   endif
 
 endfunction
@@ -115,6 +142,9 @@
 %! assert (std (x), [0, 0]);
 %! assert (std (y), sqrt (2), sqrt (eps));
 %! assert (std (x, 0, 2), zeros (10, 1));
+%! assert (std (x, [1:10]), [0, 0]);
+%! assert (std (y, [1, 4]), .8);
+%! assert (std ([1:10; 10:2:28]', [1:5, 10:14]), [2.2647, 4.5294], 1e-4);
 
 %!assert (std (ones (3, 1, 2), 0, 2), zeros (3, 1, 2))
 %!assert (std ([1 2], 0), sqrt (2)/2, 5*eps)
@@ -128,7 +158,9 @@
 ## Test input validation
 %!error <Invalid call> std ()
 %!error <X must be a numeric> std (['A'; 'B'])
-%!error <OPT must be 0 or 1> std (1, 2)
+%!error <W must be 0 or 1> std (1, 2)
+%!error <W must be .* a vector of nonnegative scalars> std ([1 2], [-1 0])
+%!error <the length of W must be> std ([1 2], [1 2 3])
 %!error <DIM must be an integer> std (1, [], ones (2,2))
 %!error <DIM must be an integer> std (1, [], 1.5)
 %!error <DIM must be .* a valid dimension> std (1, [], 0)
