# HG changeset patch
# User Stefano Guidoni <ilguido@users.sf.net>
# Date 1635002316 -7200
#      Sat Oct 23 17:18:36 2021 +0200
# Node ID 48d897408527292f3535fd3463e1f187b6926a77
# Parent  ba4aebad10d1c40d160e466361ad8826c4893361
var.m: add weighted variance (patch #10103)

diff -r ba4aebad10d1 -r 48d897408527 NEWS
--- a/NEWS	Wed Oct 13 10:36:20 2021 -0700
+++ b/NEWS	Sat Oct 23 17:18:36 2021 +0200
@@ -244,6 +244,9 @@
 compatibility the `print_struct_array_contents` internal variable is set
 to true.
 
+- The function `var` now accepts a weight vector as input and computes the
+weigthed variance.
+
 ### Alphabetical list of new functions added in Octave 7
 
 * `cospi`
diff -r ba4aebad10d1 -r 48d897408527 scripts/statistics/var.m
--- a/scripts/statistics/var.m	Wed Oct 13 10:36:20 2021 -0700
+++ b/scripts/statistics/var.m	Sat Oct 23 17:18:36 2021 +0200
@@ -25,8 +25,8 @@
 
 ## -*- texinfo -*-
 ## @deftypefn  {} {} var (@var{x})
-## @deftypefnx {} {} var (@var{x}, @var{opt})
-## @deftypefnx {} {} var (@var{x}, @var{opt}, @var{dim})
+## @deftypefnx {} {} var (@var{x}, @var{w})
+## @deftypefnx {} {} var (@var{x}, @var{w}, @var{dim})
 ## Compute the variance of the elements of the vector @var{x}.
 ##
 ## The variance is defined as
@@ -53,16 +53,22 @@
 ## If @var{x} is a matrix, compute the variance 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
 ## @item 0:
-##   normalize with @math{N-1}, provides the best unbiased estimator of the
-## variance [default]
+##   normalize with @math{N-1}, provides the square root of the best unbiased
+## estimator of the variance [default]
 ##
 ## @item 1:
-##   normalizes with @math{N}, this provides the second moment around the mean
+##   normalize with @math{N}, this provides the square root of the second
+## moment around the mean
+##
+## @item a vector:
+##   compute the weighted variance 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 @math{N} is equal to 1 the value of @var{opt} is ignored and
@@ -72,7 +78,7 @@
 ## @seealso{cov, std, skewness, kurtosis, moment}
 ## @end deftypefn
 
-function retval = var (x, opt = 0, dim)
+function retval = var (x, w = 0, dim)
 
   if (nargin < 1)
     print_usage ();
@@ -82,10 +88,13 @@
     error ("var: X must be a numeric vector or matrix");
   endif
 
-  if (isempty (opt))
-    opt = 0;
-  elseif (opt != 0 && opt != 1)
-    error ("var: 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 ("var: W must be 0 or 1 or a vector of nonnegative scalars");
   endif
 
   nd = ndims (x);
@@ -99,36 +108,62 @@
     endif
   endif
 
-  n = size (x, dim);
-  if (n == 1)
-    if (isa (x, "single"))
-      retval = zeros (sz, "single");
+  if (isempty (x))
+    ## MATLAB compatibility
+    retval = NaN;
+  else
+    n = size (x, dim);
+    if (n == 1)
+      if (isa (x, "single"))
+        retval = zeros (sz, "single");
+      else
+        retval = zeros (sz);
+      endif
     else
-      retval = zeros (sz);
+      if (isscalar (w))
+        retval = sumsq (center (x, dim), dim) / (n - 1 + w);
+      else
+        ## Weighted variance
+        if (length (w) != n)
+          error (["var: the length of W must be equal to the size of X "...
+                  "along which the variance 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, dim) ./ sum (w);
+          retval = sum (w .* ((x .- mu) .^ 2), dim) / den;
+        endif
+      endif
     endif
-  elseif (numel (x) > 0)
-    retval = sumsq (center (x, dim), dim) / (n - 1 + opt);
-  else
-    error ("var: X must not be empty");
   endif
 
 endfunction
 
 
+%!assert (var ([]), NaN)
 %!assert (var (13), 0)
 %!assert (var (single (13)), single (0))
 %!assert (var ([1,2,3]), 1)
 %!assert (var ([1,2,3], 1), 2/3, eps)
 %!assert (var ([1,2,3], [], 1), [0,0,0])
 %!assert (var ([1,2,3], [], 3), [0,0,0])
+%!assert (var ([1:7], [1:7]), 3)
+%!assert (var ([(eye (3))], [1:3]), [.1389 .2222 .25], 1e-4)
+%!assert (var (ones (2,2,2), [1:2], 3), [(zeros (2,2))])
 
 ## Test input validation
 %!error <Invalid call> var ()
 %!error <X must be a numeric> var (['A'; 'B'])
-%!error <OPT must be 0 or 1> var (1, -1)
-%!error <FLAG must be 0 or 1> skewness (1, 2)
-%!error <FLAG must be 0 or 1> skewness (1, [1 0])
+%!error <W must be 0 or 1> var (1, 2)
+%!error <W must be .* a vector of nonnegative scalars> var ([1 2], [-1 0])
+%!error <W must be .* a vector of nonnegative scalars> var ([1 2], eye (2))
+%!error <the length of W must be> var ([1 2], [1 2 3])
 %!error <DIM must be an integer> var (1, [], ones (2,2))
 %!error <DIM must be an integer> var (1, [], 1.5)
 %!error <DIM must be .* a valid dimension> var (1, [], 0)
-%!error <X must not be empty> var ([], 1)
