# HG changeset patch
# User Nicholas R. Jankowski <jankowski.nicholas@gmail.com>
# Date 1652375452 14400
# Node ID 5330efaf9476927ca49b5168fa2abd91e97a718b
# Parent  4298af839d2025c192bdaa726095a4f1903b86bb
Add optional second output to var and std (bug #62395)

* scripts/statistics/var.m: Add optional second output containing the mean
used to calculate the variance.  Move weight isempty check ahead of vector
dimension isscalar check to avoid triggering incompatability error.  Add BISTs
testing second output with different calling options.  Add BIST testing empty
value passed as variance weight treated as zero.  Add new output behavior to
docstring, and update function definitions to show the primary variable.
* scripts/statistics/std.m: Add passthrough for second output from var when
std called with two outputs.  Add BISTs testing second output with different
calling options. Update docstring noting new output behavior.
* etc/NEWS.8.md: Note output changes to var and std under Matlab Compatability.

diff -r 4298af839d20 -r 5330efaf9476 etc/NEWS.8.md
--- a/etc/NEWS.8.md	Wed May 11 12:14:08 2022 -0700
+++ b/etc/NEWS.8.md	Thu May 12 13:10:52 2022 -0400
@@ -36,10 +36,12 @@
 - `format` now accepts the option "default", which is equivalent to
   calling `format` without any options to reset the default state.
 
-
 - `quadgk` now stops iterating when `error <= tolerance` while the previous
   condition was `error < tolerance`.
 
+- `var` and `std` now optionally output a second argument containing the mean
+  or weighted mean.
+
 ### Alphabetical list of new functions added in Octave 8
 
 
diff -r 4298af839d20 -r 5330efaf9476 scripts/statistics/std.m
--- a/scripts/statistics/std.m	Wed May 11 12:14:08 2022 -0700
+++ b/scripts/statistics/std.m	Thu May 12 13:10:52 2022 -0400
@@ -28,6 +28,7 @@
 ## @deftypefnx {} {@var{y} =} std (@var{x}, @var{w})
 ## @deftypefnx {} {@var{y} =} std (@var{x}, @var{w}, @var{dim})
 ## @deftypefnx {} {@var{y} =} std (@var{x}, @var{w}, @qcode{"ALL"})
+## @deftypefnx {} {[@var{y}, @var{mu}] =} std (@dots{})
 ## Compute the standard deviation of the elements of the vector @var{x}.
 ##
 ## The standard deviation is defined as
@@ -84,12 +85,21 @@
 ## elements of @var{x}, and is equivalent to @code{std (@var{x}(:))}.
 ##
 ## When @var{dim} is a vector or @qcode{"ALL"}, @var{w} must be either 0 or 1.
+##
+## If requested the optional second output variable @var{mu} will contain the
+## mean or weighted mean used to calcluate @var{y}, and will be the same size
+## as @var{y}.
 ## @seealso{var, bounds, mad, range, iqr, mean, median}
 ## @end deftypefn
 
-function y = std (varargin)
+function [y, mu] = std (varargin)
 
-  y = sqrt (var (varargin{:}));
+  if (nargout < 2)
+    y = sqrt (var (varargin{:}));
+  else
+    [y, mu] = var (varargin{:});
+    y = sqrt (y);
+  endif
 
 endfunction
 
@@ -128,6 +138,27 @@
 %!assert (std (ones (1,3,0,2), [], 3), NaN(1,3,1,2))
 %!assert (std (ones (1,3,0,2), [], 4), NaN(1,3,0))
 
+##Test second output
+%!test <*62395>
+%! [~, m] = std (1);
+%! assert (m, 1);
+%! [~, m] = std ([1,2,3; 3,2,1]);
+%! assert (m, [2,2,2]);
+%! [~, m] = std ([]);
+%! assert (m, NaN);
+%! [~, m] = std ([1 2 3], 0);
+%! assert (m, 2);
+%! [~, m] = std ([1 2 3], 1);
+%! assert (m, 2);
+%! [~, m] = std ([1 2 3], [1 2 3]);
+%! assert (m, 7/3, eps);
+%! [~, m] = std ([1 2 3], [], 1);
+%! assert (m, [1 2 3]);
+%! [~, m] = std ([1 2 3; 1,2,3], [], [1 2]);
+%! assert (m, 2);
+%! [~, m] = std ([1 2 3; 1,2,3], [], "all");
+%! assert (m, 2);
+
 
 ## Test input validation
 %!error <Invalid call> std ()
diff -r 4298af839d20 -r 5330efaf9476 scripts/statistics/var.m
--- a/scripts/statistics/var.m	Wed May 11 12:14:08 2022 -0700
+++ b/scripts/statistics/var.m	Thu May 12 13:10:52 2022 -0400
@@ -24,10 +24,11 @@
 ########################################################################
 
 ## -*- texinfo -*-
-## @deftypefn  {} {} var (@var{x})
-## @deftypefnx {} {} var (@var{x}, @var{w})
-## @deftypefnx {} {} var (@var{x}, @var{w}, @var{dim})
-## @deftypefnx {} {} var (@var{x}, @var{w}, @qcode{"ALL"})
+## @deftypefn  {} {@var{v} =} var (@var{x})
+## @deftypefnx {} {@var{v} =} var (@var{x}, @var{w})
+## @deftypefnx {} {@var{v} =} var (@var{x}, @var{w}, @var{dim})
+## @deftypefnx {} {@var{v} =} var (@var{x}, @var{w}, @qcode{"ALL"})
+## @deftypefnx {} {[@var{v}, @var{m}] =} var (@dots{})
 ## Compute the variance of the elements of the vector @var{x}.
 ##
 ## The variance is defined as
@@ -85,10 +86,14 @@
 ## elements of @var{x}, and is equivalent to @code{var (@var{x}(:))}.
 ##
 ## When @var{dim} is a vector or @qcode{"ALL"}, @var{w} must be either 0 or 1.
+##
+## If requested the optional second output variable @var{mu} will contain the
+## mean or weighted mean used to calcluate @var{v}, and will be the same size
+## as @var{v}.
 ## @seealso{cov, std, skewness, kurtosis, moment}
 ## @end deftypefn
 
-function retval = var (x, w = 0, dim)
+function [v, mu] = var (x, w = 0, dim)
 
   if (nargin < 1)
     print_usage ();
@@ -100,12 +105,17 @@
     error ("var: X must be a numeric vector or matrix");
   endif
 
+  if (isempty (w))
+    w = 0;
+  endif
+
   nd = ndims (x);
   sz = size (x);
   emptydimflag = false;
 
   if (isempty (dim))
-    emptydimflag = true;  ## Compatibliity hack for empty x, ndims==2
+    emptydimflag = true;  # Compatibliity hack for empty x, ndims==2
+
     ## Find the first non-singleton dimension.
    (dim = find (sz != 1, 1)) || (dim = 1);
 
@@ -169,9 +179,7 @@
   endif
 
   n = size (x, dim);
-  if (isempty (w))
-    w = 0;
-  elseif (! isvector (w) ||
+  if (! isvector (w) ||
           ! isnumeric (w) ||
           (isvector (w) && any (w < 0)) ||
           (isscalar (w) && ((w != 0 && w != 1) && (n != 1))))
@@ -180,11 +188,13 @@
 
   if (isempty (x))
     if (emptydimflag && isequal (sz, [0 0]))
-      retval = NaN;
+      v = NaN;
+      mu = NaN;
     else
       output_size = sz;
       output_size(dim) = 1;
-      retval = NaN(output_size);
+      v = NaN(output_size);
+      mu = NaN(output_size);
     endif
   else
     if (n == 1)
@@ -193,14 +203,19 @@
                   "in the dimension along which variance is calculated"])
       else
         if (isa (x, "single"))
-          retval = zeros (sz, "single");
+          v = zeros (sz, "single");
+          mu = x;
         else
-          retval = zeros (sz);
+          v = zeros (sz);
+          mu = x;
         endif
       endif
     else
       if (isscalar (w))
-        retval = sumsq (center (x, dim), dim) / (n - 1 + w);
+        v = sumsq (center (x, dim), dim) / (n - 1 + w);
+        if (nargout == 2)
+          mu = mean (x, dim);
+        endif
       else
         ## Weighted variance
         if (length (w) != n)
@@ -215,8 +230,8 @@
             w = reshape (w, newdims);
           endif
           den = sum (w);
-          mu = sum (w .* x, dim) ./ sum (w);
-          retval = sum (w .* ((x - mu) .^ 2), dim) / den;
+          mu = sum (w .* x, dim) ./ den;
+          v = sum (w .* ((x - mu) .^ 2), dim) ./ den;
         endif
       endif
     endif
@@ -238,6 +253,7 @@
 %!assert (var (ones (2,2,2), [1:2], 3), [(zeros (2,2))])
 %!assert (var ([1 2; 3 4], 0, 'all'), var ([1:4]))
 %!assert (var (reshape ([1:8], 2, 2, 2), 0, [1 3]), [17/3 17/3], eps)
+%!assert (var ([1 2 3;1 2 3], [], [1 2]), 0.8, eps)
 
 ##Test empty inputs
 %!assert (var ([]), NaN)
@@ -259,6 +275,49 @@
 %!assert (var (ones (1,3,0,2), [], 3), NaN(1,3,1,2))
 %!assert (var (ones (1,3,0,2), [], 4), NaN(1,3,0))
 
+##Test second output
+%!test <*62395>
+%! [~, m] = var (13);
+%! assert (m, 13);
+%! [~, m] = var (single(13));
+%! assert (m, single(13));
+%! [~, m] = var ([1, 2, 3; 3 2 1], []);
+%! assert (m, [2 2 2]);
+%! [~, m] = var ([1, 2, 3; 3 2 1], [], 1);
+%! assert (m, [2 2 2]);
+%! [~, m] = var ([1, 2, 3; 3 2 1], [], 2);
+%! assert (m, [2 2]');
+%! [~, m] = var ([1, 2, 3; 3 2 1], [], 3);
+%! assert (m, [1 2 3; 3 2 1]);
+
+## 2nd output, weighted inputs, vector dims
+%!test <*62395>
+%! [~, m] = var(5,99);
+%! assert (m, 5);
+%! [~, m] = var ([1:7], [1:7]);
+%! assert (m, 5);
+%! [~, m] = var ([eye(3)], [1:3]);
+%! assert (m, [1/6, 1/3, 0.5], eps);
+%! [~, m] = var (ones (2,2,2), [1:2], 3);
+%! assert (m, ones (2,2));
+%! [~, m] = var ([1 2; 3 4], 0, 'all');
+%! assert (m, 2.5, eps);
+%! [~, m] = var (reshape ([1:8], 2, 2, 2), 0, [1 3]);
+%! assert (m, [3.5, 5.5], eps);
+
+## 2nd output, empty inputs
+%!test <*62395>
+%! [~, m] = var ([]);
+%! assert (m, NaN);
+%! [~, m] = var ([],[],1);
+%! assert (m, NaN(1,0));
+%! [~, m] = var ([],[],2);
+%! assert (m, NaN(0,1));
+%! [~, m] = var ([],[],3);
+%! assert (m, []);
+%! [~, m] = var (ones (1,3,0,2));
+%! assert (m, NaN(1,1,0,2));
+
 ## Test input validation
 %!error <Invalid call> var ()
 %!error <X must be a numeric> var (['A'; 'B'])