Mercurial > octave
changeset 31868:c6eeb8b44c28
Port statistics pkg mean, median, var, & std into core (patch #10314)
* mean.m: Overhaul function to vectorize 'nanflag' codepath, improve 'vecdim'
handling for arbitrary dimensions (bug #63460), ensure proper outputs for
empty inputs (bug #50583), and avoid overflow when summing large int values or
precision loss for trying to handle large int64/uint64 values as doubles
(bug #54567). Add BISTs for items above and to address corner case inputs.
* median.m: Overhaul function to implement 'nanflag' option to optionally
ignore NaN values (bug #50571), implement 'all' (bug #58116) and 'VECDIM'
(bug #58089) dimension input options, ensure proper outputs for empty inputs
(bug #50583), and avoid overflow when summing large int values and precision
loss for trying to handle large int64/uint64 values as doubles (bug #54567).
Add BISTs for items above and to address corner case inputs.
* var.m: Overhaul function to implement 'nanflag' option to optionally ignore
NaN values (bug #50571), implement 'all' (bug #58116) and 'VECDIM' (bug #58089)
dimension input options, avoid broadcasting operations from erroring on sparse
or diagonal type inputs (bug #63291), ensure proper outputs for empty inputs
(bug #50583). Add BISTs for items above and to address corner case inputs.
* std.m: Update docstring to match changes in var.
* NEWS.9.md: Note changes to above functions under Matlab Compatibility.
| author | Nicholas R. Jankowski <jankowski.nicholas@gmail.com> |
|---|---|
| date | Wed, 01 Mar 2023 00:38:54 -0500 |
| parents | a9c8b1f8fb32 |
| children | 0149a2f4ac98 |
| files | etc/NEWS.9.md scripts/statistics/mean.m scripts/statistics/median.m scripts/statistics/std.m scripts/statistics/var.m |
| diffstat | 5 files changed, 1774 insertions(+), 414 deletions(-) [+] |
line wrap: on
line diff
--- a/etc/NEWS.9.md Tue Feb 28 15:45:03 2023 -0500 +++ b/etc/NEWS.9.md Wed Mar 01 00:38:54 2023 -0500 @@ -29,6 +29,14 @@ ### Matlab compatibility +- Overhauled `mean`, `median`, `var`, and `std` functions have been imported +from statistics package v1.5.4 to compatibly implement 'nanflag' (bug #50571), +'all' (bug #58116), and 'vecdim' (bug #58089) options, preserve output class, +and match expected output behavior for empty (bug #50583) and sparse inputs +(bug #63291). Both `median` and `mean` can handle large int values without +overflow or precision concerns (bug #54567). `median` also adopts the +'outtype' option from `mean`. + ### Alphabetical list of new functions added in Octave 9 * `tensorprod`
--- a/scripts/statistics/mean.m Tue Feb 28 15:45:03 2023 -0500 +++ b/scripts/statistics/mean.m Wed Mar 01 00:38:54 2023 -0500 @@ -1,6 +1,6 @@ ######################################################################## ## -## Copyright (C) 1995-2023 The Octave Project Developers +## Copyright (C) 1996-2023 The Octave Project Developers ## ## See the file COPYRIGHT.md in the top-level directory of this ## distribution or <https://octave.org/copyright/>. @@ -24,21 +24,19 @@ ######################################################################## ## -*- texinfo -*- -## @deftypefn {} {@var{y} =} mean (@var{x}) -## @deftypefnx {} {@var{y} =} mean (@var{x}, 'all') -## @deftypefnx {} {@var{y} =} mean (@var{x}, @var{dim}) -## @deftypefnx {} {@var{y} =} mean (@dots{}, '@var{outtype}') -## @deftypefnx {} {@var{y} =} mean (@dots{}, '@var{nanflag}') +## @deftypefn {} {@var{m} =} mean (@var{x}) +## @deftypefnx {} {@var{m} =} mean (@var{x}, @var{dim}) +## @deftypefnx {} {@var{m} =} mean (@var{x}, @var{vecdim}) +## @deftypefnx {} {@var{m} =} mean (@var{x}, "all") +## @deftypefnx {} {@var{m} =} mean (@dots{}, @var{nanflag}) +## @deftypefnx {} {@var{m} =} mean (@dots{}, @var{outtype}) ## Compute the mean of the elements of @var{x}. ## -## @itemize -## @item ## If @var{x} is a vector, then @code{mean (@var{x})} returns the ## mean of the elements in @var{x} defined as ## @tex ## $$ {\rm mean}(x) = \bar{x} = {1\over N} \sum_{i=1}^N x_i $$ ## where $N$ is the number of elements of @var{x}. -## ## @end tex ## @ifnottex ## @@ -47,29 +45,29 @@ ## @end example ## ## @noindent -## where @math{N} is the number of elements in the @var{x} vector. +## where @math{N} is the number of elements in @var{x}. ## ## @end ifnottex ## -## @item -## If @var{x} is a matrix, then @code{mean} returns a row vector with the mean -## of each column in @var{x}. +## If @var{x} is an array, then @code{mean(@var{x})} computes the mean along the +## first nonsingleton dimension of @var{x}. +## +## The optional variable @var{dim} forces @code{mean} to operate over the +## specified dimension, which must be a positive integer-valued number. +## Specifying any singleton dimension in @var{x}, including any dimension +## exceeding @code{ndims (@var{x})}, will result in a mean equal to @var{x}. ## -## @item -## If @var{x} is a multi-dimensional array, then @code{mean} operates along the -## first non-singleton dimension of @var{x}. -## @end itemize +## Specifying the dimensions as @var{vecdim}, a vector of non-repeating +## dimensions, will return the mean over the array slice defined by +## @var{vecdim}. If @var{vecdim} indexes all dimensions of @var{x}, then it is +## equivalent to the option @qcode{"all"}. Any dimension in @var{vecdim} greater +## than @code{ndims (@var{x})} is ignored. ## -## The optional input @var{dim} forces @code{mean} to operate over the -## specified dimension(s). @var{dim} can either be a scalar dimension or a -## vector of non-repeating dimensions. Dimensions must be positive integers, -## and the mean is calculated over the array slice defined by @var{dim}. -## -## Specifying dimension @qcode{"all"} will force @code{mean} to operate on all -## elements of @var{x}, and is equivalent to @code{mean (@var{x}(:))}. +## Specifying the dimension as @qcode{"all"} will force @code{mean} to operate +## on all elements of @var{x}, and is equivalent to @code{mean (@var{x}(:))}. ## ## The optional input @var{outtype} specifies the data type that is returned. -## Valid values are: +## @var{outtype} can take the following values: ## ## @table @asis ## @item @qcode{'default'} : Output is of type double, unless the input is @@ -77,174 +75,306 @@ ## ## @item @qcode{'double'} : Output is of type double. ## -## @item @qcode{'native'} : Output is of the same type as the input -## (@code{class (@var{x})}), unless the input is logical in which case the +## @item @qcode{'native'} : Output is of the same type as the input as reported +## by (@code{class (@var{x})}), unless the input is logical in which case the ## output is of type double. -## ## @end table ## -## The optional input @var{nanflag} specifies whether to include/exclude NaN -## values from the calculation. By default, NaN values are included in the -## calculation (@var{nanflag} has the value @qcode{'includenan'}). To exclude -## NaN values, set the value of @var{nanflag} to @qcode{'omitnan'}. +## The optional variable @var{nanflag} specifies whether to include or exclude +## NaN values from the calculation using any of the previously specified input +## argument combinations. The default value for @var{nanflag} is "includenan" +## which keeps NaN values in the calculation. To exclude NaN values set the +## value of @var{nanflag} to "omitnan". The output will still contain NaN +## values if @var{x} consists of all NaN values in the operating dimension. ## -## @seealso{median, mode} +## @seealso{median, mode, movmean} ## @end deftypefn -function y = mean (x, varargin) +function m = mean (x, varargin) + + if (nargin < 1 || nargin > 4) + print_usage (); + endif - if (nargin < 1 || nargin > 4 || ! all (cellfun (@ischar, varargin(2:end)))) + ## Set initial conditions + all_flag = 0; + omitnan = 0; + out_flag = 0; + + nvarg = numel (varargin); + varg_chars = cellfun ("ischar", varargin); + outtype = "default"; + szx = size (x); + ndx = ndims (x); + + if (nvarg > 1 && ! varg_chars(2:end)) + ## Only first varargin can be numeric print_usage (); endif - ## Check all char arguments. - all_flag = false; - omitnan = false; - outtype = "default"; - - for i = 1:numel (varargin) - if (ischar (varargin{i})) - switch (varargin{i}) + ## Process any other char arguments. + if (any (varg_chars)) + for i = varargin(varg_chars) + switch (lower (i{:})) case "all" all_flag = true; + + case "omitnan" + omitnan = true; + case "includenan" omitnan = false; - case "omitnan" - omitnan = true; - case {"default", "double", "native"} - outtype = varargin{i}; + + case "default" + if (out_flag) + error ("mean: only one OUTTYPE can be specified.") + endif + if (isa (x, "single")) + outtype = "single"; + else + outtype = "double"; + endif + out_flag = 1; + + case "native" + outtype = class (x); + if (out_flag) + error ("mean: only one OUTTYPE can be specified.") + elseif (strcmp (outtype, "logical")) + outtype = "double"; + elseif (strcmp (outtype, "char")) + error ("mean: OUTTYPE 'native' cannot be used with char inputs."); + endif + out_flag = 1; + + case "double" + if (out_flag) + error ("mean: only one OUTTYPE can be specified.") + endif + outtype = "double"; + out_flag = 1; + otherwise print_usage (); endswitch + endfor + varargin(varg_chars) = []; + nvarg = numel (varargin); + endif + + if strcmp (outtype, "default") + if (isa (x, "single")) + outtype = "single"; + else + outtype = "double"; endif - endfor - varargin(cellfun (@ischar, varargin)) = []; + endif - if (((numel (varargin) == 1) && ! (isnumeric (varargin{1}))) ... - || (numel (varargin) > 1)) + if ((nvarg > 1) || ((nvarg == 1) && ! (isnumeric (varargin{1})))) + ## After trimming char inputs can only be one varargin left, must be numeric print_usage (); endif - if (! (isnumeric (x) || islogical (x))) - error ("mean: X must be either a numeric or logical vector or matrix"); + if (! (isnumeric (x) || islogical (x) || ischar (x))) + error ("mean: X must be either a numeric, boolean, or character array."); endif - if (numel (varargin) == 0) - + ## Process special cases for in/out size + if (nvarg == 0) ## Single numeric input argument, no dimensions given. + if (all_flag) - n = numel (x(:)); + x = x(:); + if (omitnan) - idx = isnan (x); - n -= sum (idx(:)); - x(idx) = 0; + x = x(! isnan (x)); endif - y = sum (x(:), 1) ./ n; + + if (any (isa (x, {"int64", "uint64"}))) + m = int64_mean (x, 1, numel (x), outtype); + else + m = sum (x) ./ numel (x); + endif + else - sz = size (x); ## Find the first non-singleton dimension. - (dim = find (sz != 1, 1)) || (dim = 1); - n = size (x, dim); + (dim = find (szx != 1, 1)) || (dim = 1); + n = szx(dim); if (omitnan) idx = isnan (x); n = sum (! idx, dim); x(idx) = 0; endif - y = sum (x, dim) ./ n; + + if (any (isa (x, {"int64", "uint64"}))) + m = int64_mean (x, dim, n, outtype); + else + m = sum (x, dim) ./ n; + endif + endif else ## Two numeric input arguments, dimensions given. Note scalar is vector! - dim = varargin{1}; - if (! (isvector (dim) && all (dim > 0) && all (rem (dim, 1) == 0))) - error ("mean: DIM must be a positive integer scalar or vector"); + vecdim = varargin{1}; + if (isempty (vecdim) || ! (isvector (vecdim) && all (vecdim > 0)) ... + || any (rem (vecdim, 1))) + error ("mean: DIM must be a positive integer scalar or vector."); endif - if (isscalar (dim)) - - n = size (x, dim); - if (omitnan) - idx = isnan (x); - n = sum (! idx, dim); - x(idx) = 0; - endif - y = sum (x, dim) ./ n; - + if (ndx == 2 && isempty (x) && szx == [0,0]) + ## FIXME: this special case handling could be removed once sum + ## compatably handles all sizes of empty inputs + sz_out = szx; + sz_out (vecdim(vecdim <= ndx)) = 1; + m = NaN (sz_out); else - sz = size (x); - ndims = numel (sz); - misdim = [1:ndims]; - - dim(dim > ndims) = []; # weed out dimensions larger than array - misdim(dim) = []; # remove dims asked for leaving missing dims - - switch (numel (misdim)) - ## if all dimensions are given, compute x(:) - case 0 - n = numel (x(:)); + if (isscalar (vecdim)) + if (vecdim > ndx) + m = x; + else + n = szx(vecdim); if (omitnan) - idx = isnan (x); - n -= sum (idx(:)); - x(idx) = 0; + nanx = isnan (x); + n = sum (! nanx, vecdim); + x(nanx) = 0; endif - y = sum (x(:), 1) ./ n; + + if (any (isa (x, {"int64", "uint64"}))) + m = int64_mean (x, vecdim, n, outtype); + else + m = sum (x, vecdim) ./ n; + endif + + endif - ## for 1 dimension left, return column vector - case 1 - x = permute (x, [misdim, dim]); - y = zeros (size (x, 1), 1, "like", x); - for i = 1:size (x, 1) - x_vec = x(i,:)(:); + else + vecdim = sort (vecdim); + if (! all (diff (vecdim))) + error ("mean: VECDIM must contain non-repeating positive integers."); + endif + ## Ignore exceeding dimensions in VECDIM + vecdim(find (vecdim > ndims (x))) = []; + + if (isempty (vecdim)) + m = x; + else + ## Move vecdims to dim 1. + + ## Calculate permutation vector + remdims = 1 : ndx; # All dimensions + remdims(vecdim) = []; # Delete dimensions specified by vecdim + nremd = numel (remdims); + + ## If all dimensions are given, it is similar to all flag + if (nremd == 0) + x = x(:); if (omitnan) - x_vec = x_vec(! isnan (x_vec)); + x = x(! isnan (x)); + endif + + if (any (isa (x, {"int64", "uint64"}))) + m = int64_mean (x, 1, numel (x), outtype); + else + m = sum (x) ./ numel (x); endif - y(i) = sum (x_vec, 1) ./ numel (x_vec); - endfor - y = ipermute (y, [misdim, dim]); + + else + ## Permute to bring vecdims to front + perm = [vecdim, remdims]; + x = permute (x, perm); + + ## Reshape to squash all vecdims in dim1 + num_dim = prod (szx(vecdim)); + szx(vecdim) = []; + szx = [ones(1, numel(vecdim)), szx]; + szx(1) = num_dim; + x = reshape (x, szx); - ## for 2 dimensions left, return matrix - case 2 - x = permute (x, [misdim, dim]); - y = zeros (size (x, 1), size (x, 2), "like", x); - for i = 1:size (x, 1) - for j = 1:size (x, 2) - x_vec = x(i,j,:)(:); - if (omitnan) - x_vec = x_vec(! isnan (x_vec)); - endif - y(i,j) = sum (x_vec, 1) ./ numel (x_vec); - endfor - endfor - y = ipermute (y, [misdim, dim]); + ## Calculate mean on dim1 + if (omitnan) + nanx = isnan (x); + x(nanx) = 0; + n = sum (! nanx, 1); + else + n = szx(1); + endif - ## for more than 2 dimensions left, throw error - otherwise - error ("DIM must index at least N-2 dimensions of X"); - endswitch + if (any (isa (x, {"int64", "uint64"}))) + m = int64_mean (x, 1, n, outtype); + else + m = sum (x, 1) ./ n; + endif + + ## Inverse permute back to correct dimensions + m = ipermute (m, perm); + endif + endif + endif endif - endif - ## Convert output if requested - switch (outtype) - case "default" - ## do nothing, the operators already do the right thing. - case "double" - y = double (y); - case "native" - if (! islogical (x)) - y = cast (y, class (x)); - endif - otherwise - ## FIXME: This is unreachable code. Valid values already - ## checked in input validation. - error ("mean: OUTTYPE '%s' not recognized", outtype); - endswitch - + ## Convert output as requested + if (! strcmp (class (m), outtype)) + switch (outtype) + case "double" + m = double (m); + case "single" + m = single (m); + otherwise + if (! islogical (x)) + m = cast (m, outtype); + endif + endswitch + endif endfunction +function m = int64_mean (x, dim, n, outtype) + ## Avoid int overflow in large ints. Smaller ints processed as double + ## avoids overflow but large int64 values have floating pt error as double. + ## Use integer math and manual remainder correction to avoid this. + if (any (abs (x(:)) >= flintmax / n)) + rmdr = double (rem (x, n)) / n; + rmdr_hilo = logical (int8 (rmdr)); # Integer rounding direction indicator + + ## Do 'native' int summation to prevent double precision error, + ## then add back in lost round-up/down remainders. + + m = sum (x/n, dim, "native"); + + ## rmdr.*!rmdr_hilo = remainders that were rounded down in abs val + ## signs retained, can be summed and added back. + ## rmdr.*rmdr_hilo = remainders that were rounded up in abs val. + ## need to add back difference between 1 and rmdr, retaining sign. + + rmdr = sum (rmdr .* !rmdr_hilo, dim) - ... + sum ((1 - abs (rmdr)) .* rmdr_hilo .* sign (rmdr), dim); + + if (any (abs (m(:)) >= flintmax)) + ## Avoid float errors when combining for large m. + ## FIXME - may also need to include checking rmdir for large numel (x), + ## as its value could be on the order of numel (x). + if (any (strcmp (outtype, {"int64", "uint64"}))) + m = m + rmdr; + else + m = double (m) + rmdr; + endif + + else + m = double(m) + rmdr; + switch (outtype) + case "int64" + m = int64 (m); + case "uint64" + m = uint64 (m); + endswitch + endif + else + m = double (sum (x, dim, "native")) ./ n; + endif +endfunction %!test %! x = -10:10; @@ -260,34 +390,121 @@ %!assert (mean (single ([1 0 1 1])), single (0.75)) %!assert (mean ([1 2], 3), [1 2]) -## Test outtype option +#### Test outtype option %!test %! in = [1 2 3]; %! out = 2; %! assert (mean (in, "default"), mean (in)); %! assert (mean (in, "default"), out); -%! +%! assert (mean (in, "double"), out); +%! assert (mean (in, "native"), out); + +%!test %! in = single ([1 2 3]); %! out = 2; %! assert (mean (in, "default"), mean (in)); %! assert (mean (in, "default"), single (out)); %! assert (mean (in, "double"), out); %! assert (mean (in, "native"), single (out)); -%! + +%!test +%! in = logical ([1 0 1]); +%! out = 2/3; +%! assert (mean (in, "default"), mean (in), eps); +%! assert (mean (in, "default"), out, eps); +%! assert (mean (in, "double"), out, eps); +%! assert (mean (in, "native"), out, eps); + +%!test +%! in = char ("ab"); +%! out = 97.5; +%! assert (mean (in, "default"), mean (in), eps); +%! assert (mean (in, "default"), out, eps); +%! assert (mean (in, "double"), out, eps); + +%!test %! in = uint8 ([1 2 3]); %! out = 2; %! assert (mean (in, "default"), mean (in)); %! assert (mean (in, "default"), out); %! assert (mean (in, "double"), out); %! assert (mean (in, "native"), uint8 (out)); -%! -%! in = logical ([1 0 1]); -%! out = 2/3; + +%!test +%! in = uint8 ([0 1 2 3]); +%! out = 1.5; +%! out_u8 = 2; +%! assert (mean (in, "default"), mean (in), eps); +%! assert (mean (in, "default"), out, eps); +%! assert (mean (in, "double"), out, eps); +%! assert (mean (in, "native"), uint8 (out_u8)); +%! assert (class (mean (in, "native")), "uint8"); + +%!test ## internal sum exceeding intmax +%! in = uint8 ([3 141 141 255]); +%! out = 135; +%! assert (mean (in, "default"), mean (in)); +%! assert (mean (in, "default"), out); +%! assert (mean (in, "double"), out); +%! assert (mean (in, "native"), uint8 (out)); +%! assert (class (mean (in, "native")), "uint8"); + +%!test ## fractional answer with interal sum exceeding intmax +%! in = uint8 ([1 141 141 255]); +%! out = 134.5; +%! out_u8 = 135; %! assert (mean (in, "default"), mean (in)); %! assert (mean (in, "default"), out); -%! assert (mean (in, "native"), out); # logical ignores native option +%! assert (mean (in, "double"), out); +%! assert (mean (in, "native"), uint8 (out_u8)); +%! assert (class (mean (in, "native")), "uint8"); -## Test single input and optional arguments "all", DIM, "omitnan") +%!test <54567> ## large int64 sum exceeding intmax and double precision limit +%! in_same = uint64 ([intmax("uint64") intmax("uint64")-2]); +%! out_same = intmax ("uint64")-1; +%! in_opp = int64 ([intmin("int64"), intmax("int64")-1]); +%! out_opp = -1; +%! in_neg = int64 ([intmin("int64") intmin("int64")+2]); +%! out_neg = intmin ("int64")+1; +%! +%! ## both positive +%! assert (mean (in_same, "default"), mean (in_same)); +%! assert (mean (in_same, "default"), double (out_same)); +%! assert (mean (in_same, "double"), double (out_same)); +%! assert (mean (in_same, "native"), uint64 (out_same)); +%! assert (class (mean (in_same, "native")), "uint64"); +%! +%! ## opposite signs +%! assert (mean (in_opp, "default"), mean (in_opp)); +%! assert (mean (in_opp, "default"), double (out_opp)); +%! assert (mean (in_opp, "double"), double (out_opp)); +%! assert (mean (in_opp, "native"), int64 (out_opp)); +%! assert (class (mean (in_opp, "native")), "int64"); +%! +%! ## both negative +%! assert (mean (in_neg, "default"), mean (in_neg)); +%! assert (mean (in_neg, "default"), double(out_neg)); +%! assert (mean (in_neg, "double"), double(out_neg)); +%! assert (mean (in_neg, "native"), int64(out_neg)); +%! assert (class (mean (in_neg, "native")), "int64"); + +## Additional tests int64 and double precision limits +%!test <54567> +%! in = [(intmin('int64')+5), (intmax('int64'))-5]; +%! assert (mean (in, "native"), int64(-1)); +%! assert (class (mean (in, "native")), "int64"); +%! assert (mean (double(in)), double(0) ); +%! assert (mean (in), double(-0.5) ); +%! assert (mean (in, "default"), double(-0.5) ); +%! assert (mean (in, "double"), double(-0.5) ); +%! assert (mean (in, "all", "native"), int64(-1)); +%! assert (mean (in, 2, "native"), int64(-1)); +%! assert (mean (in, [1 2], "native"), int64(-1)); +%! assert (mean (in, [2 3], "native"), int64(-1)); +%! assert (mean ([intmin("int64"), in, intmax("int64")]), double(-0.5)) +%! assert (mean ([in; int64([1 3])], 2, "native"), int64([-1; 2])); + +## Test input and optional arguments "all", DIM, "omitnan") %!test %! x = [-10:10]; %! y = [x;x+5;x-5]; @@ -298,6 +515,8 @@ %! assert (mean (y', "omitnan"), [0 5.35 -5]); %! z = y + 20; %! assert (mean (z, "all"), NaN); +%! assert (mean (z, "all", "includenan"), NaN); +%! assert (mean (z, "all", "omitnan"), 20.03225806451613, 4e-14); %! m = [20 NaN 15]; %! assert (mean (z'), m); %! assert (mean (z', "includenan"), m); @@ -307,7 +526,7 @@ %! assert (mean (z, 2, "native", "omitnan"), m'); %! assert (mean (z, 2, "omitnan", "native"), m'); -## Test boolean input +# Test boolean input %!test %! assert (mean (true, "all"), 1); %! assert (mean (false), 0); @@ -318,7 +537,49 @@ %! assert (mean ([true false NaN], 2, "omitnan"), 0.5); %! assert (mean ([true false NaN], 2, "omitnan", "native"), 0.5); -## Test dimension indexing with vecdim in N-dimensional arrays +## Test char inputs +%!assert (mean ("abc"), double (98)) +%!assert (mean ("ab"), double (97.5), eps) +%!assert (mean ("abc", "double"), double (98)) +%!assert (mean ("abc", "default"), double (98)) + +## Test NaN inputs +%!test +%! x = magic (4); +%! x([2, 9:12]) = NaN; +%! assert (mean (x), [NaN 8.5, NaN, 8.5], eps); +%! assert (mean (x,1), [NaN 8.5, NaN, 8.5], eps); +%! assert (mean (x,2), NaN(4,1), eps); +%! assert (mean (x,3), x, eps); +%! assert (mean (x, 'omitnan'), [29/3, 8.5, NaN, 8.5], eps); +%! assert (mean (x, 1, 'omitnan'), [29/3, 8.5, NaN, 8.5], eps); +%! assert (mean (x, 2, 'omitnan'), [31/3; 9.5; 28/3; 19/3], eps); +%! assert (mean (x, 3, 'omitnan'), x, eps); + +## Test empty inputs +%!assert (mean ([]), NaN(1,1)) +%!assert (mean (single([])), NaN(1,1,"single")) +%!assert (mean ([], 1), NaN(1,0)) +%!assert (mean ([], 2), NaN(0,1)) +%!assert (mean ([], 3), NaN(0,0)) +%!assert (mean (ones(1,0)), NaN(1,1)) +%!assert (mean (ones(1,0), 1), NaN(1,0)) +%!assert (mean (ones(1,0), 2), NaN(1,1)) +%!assert (mean (ones(1,0), 3), NaN(1,0)) +%!assert (mean (ones(0,1)), NaN(1,1)) +%!assert (mean (ones(0,1), 1), NaN(1,1)) +%!assert (mean (ones(0,1), 2), NaN(0,1)) +%!assert (mean (ones(0,1), 3), NaN(0,1)) +%!assert (mean (ones(0,1,0)), NaN(1,1,0)) +%!assert (mean (ones(0,1,0), 1), NaN(1,1,0)) +%!assert (mean (ones(0,1,0), 2), NaN(0,1,0)) +%!assert (mean (ones(0,1,0), 3), NaN(0,1,1)) +%!assert (mean (ones(0,0,1,0)), NaN(1,0,1,0)) +%!assert (mean (ones(0,0,1,0), 1), NaN(1,0,1,0)) +%!assert (mean (ones(0,0,1,0), 2), NaN(0,1,1,0)) +%!assert (mean (ones(0,0,1,0), 3), NaN(0,0,1,0)) + +## Test dimension indexing with vecdim in n-dimensional arrays %!test %! x = repmat ([1:20;6:25], [5 2 6 3]); %! assert (size (mean (x, [3 2])), [10 1 1 3]); @@ -327,33 +588,43 @@ %! assert (size (mean (x, [1 4 3])), [1 40]); %! assert (size (mean (x, [1 2 3 4])), [1 1]); -## Test results with vecdim in N-dimensional arrays and "omitnan" +## Test exceeding dimensions +%!assert (mean (ones (2,2), 3), ones (2,2)); +%!assert (mean (ones (2,2,2), 99), ones (2,2,2)); +%!assert (mean (magic (3), 3), magic (3)); +%!assert (mean (magic (3), [1 3]), [5, 5, 5]); +%!assert (mean (magic (3), [1 99]), [5, 5, 5]); + +## Test results with vecdim in n-dimensional arrays and "omitnan" %!test %! x = repmat ([1:20;6:25], [5 2 6 3]); %! m = repmat ([10.5;15.5], [5 1 1 3]); %! assert (mean (x, [3 2]), m, 4e-14); %! x(2,5,6,3) = NaN; -%! m(2,3) = NaN; +%! m(2,1,1,3) = NaN; %! assert (mean (x, [3 2]), m, 4e-14); -%! m(2,3) = 15.52301255230125; +%! m(2,1,1,3) = 15.52301255230125; %! assert (mean (x, [3 2], "omitnan"), m, 4e-14); +## Test input case insensitivity +%!assert (mean ([1 2 3], "aLL"), 2); +%!assert (mean ([1 2 3], "OmitNan"), 2); +%!assert (mean ([1 2 3], "DOUBle"), 2); + ## Test input validation -%!error <Invalid call> mean () -%!error <Invalid call> mean (1, 2, 3) -%!error <Invalid call> mean (1, 2, 3, 4, 5) -%!error <Invalid call> mean (1, "all", 3) -%!error <Invalid call> mean (1, "b") -%!error <Invalid call> mean (1, 1, "foo") -%!error <X must be either a numeric or logical> mean ({1:5}) -%!error <X must be either a numeric or logical> mean ("char") +%!error <Invalid call to mean. Correct usage is> mean () +%!error <Invalid call to mean. Correct usage is> mean (1, 2, 3) +%!error <Invalid call to mean. Correct usage is> mean (1, 2, 3, 4) +%!error <Invalid call to mean. Correct usage is> mean (1, "all", 3) +%!error <Invalid call to mean. Correct usage is> mean (1, "b") +%!error <Invalid call to mean. Correct usage is> mean (1, 1, "foo") +%!error <OUTTYPE 'native' cannot be used with char> mean ("abc", "native") +%!error <X must be either a numeric, boolean, or character> mean ({1:5}) %!error <DIM must be a positive integer> mean (1, ones (2,2)) %!error <DIM must be a positive integer> mean (1, 1.5) -%!error <DIM must be a positive integer> mean (1, -1) -%!error <DIM must be a positive integer> mean (1, -1.5) %!error <DIM must be a positive integer> mean (1, 0) -%!error <DIM must be a positive integer> mean (1, NaN) -%!error <DIM must be a positive integer> mean (1, Inf) -%!error <DIM must index at least N-2 dimensions of X> -%! mean (repmat ([1:20;6:25], [5 2 6 3 5]), [1 2]) - +%!error <DIM must be a positive integer> mean (1, []) +%!error <DIM must be a positive integer> mean (repmat ([1:20;6:25], [5 2]), -1) +%!error <DIM must be a positive integer> mean (repmat ([1:5;5:9], [5 2]), [1 -1]) +%!error <DIM must be a positive integer> mean (1, ones(1,0)) +%!error <VECDIM must contain non-repeating> mean (1, [2 2])
--- a/scripts/statistics/median.m Tue Feb 28 15:45:03 2023 -0500 +++ b/scripts/statistics/median.m Wed Mar 01 00:38:54 2023 -0500 @@ -24,9 +24,13 @@ ######################################################################## ## -*- texinfo -*- -## @deftypefn {} {@var{y} =} median (@var{x}) -## @deftypefnx {} {@var{y} =} median (@var{x}, @var{dim}) -## Compute the median value of the elements of the vector @var{x}. +## @deftypefn {} {@var{m} =} median (@var{x}) +## @deftypefnx {} {@var{m} =} median (@var{x}, @var{dim}) +## @deftypefnx {} {@var{m} =} median (@var{x}, @var{vecdim}) +## @deftypefnx {} {@var{m} =} median (@var{x}, "all") +## @deftypefnx {} {@var{m} =} median (@dots{}, @var{nanflag}) +## @deftypefnx {} {@var{m} =} median (@dots{}, @var{outtype}) +## Compute the median value of the elements of @var{x}. ## ## When the elements of @var{x} are sorted, say ## @code{@var{s} = sort (@var{x})}, the median is defined as @@ -43,65 +47,414 @@ ## ## @example ## @group -## | @var{s}(ceil(N/2)) N odd +## | @var{s}(ceil (N/2)) N odd ## median (@var{x}) = | ## | (@var{s}(N/2) + @var{s}(N/2+1))/2 N even ## @end group ## @end example ## ## @end ifnottex -## If @var{x} is of a discrete type such as integer or logical, then -## the case of even @math{N} rounds up (or toward @code{true}). +## +## If @var{x} is a array, then @code{median (@var{x})} operates along the first +## non-singleton dimension of @var{x}. +## +## The optional variable @var{dim} forces @code{median} to operate over the +## specified dimension, which must be a positive integer-valued number. +## Specifying any singleton dimension in @var{x}, including any dimension +## exceeding @code{ndims (@var{x})}, will result in a median equal to @var{x}. +## +## Specifying the dimensions as @var{vecdim}, a vector of non-repeating +## dimensions, will return the median over the array slice defined by +## @var{vecdim}. If @var{vecdim} indexes all dimensions of @var{x}, then it is +## equivalent to the option @qcode{"all"}. Any dimension in @var{vecdim} greater +## than @code{ndims (@var{x})} is ignored. +## +## Specifying the dimension as @qcode{"all"} will force @code{median} to operate +## on all elements of @var{x}, and is equivalent to @code{median (@var{x}(:))}. +## +## @code{median (@dots{}, @var{outtype})} returns the median with a specified +## data type, using any of the input arguments in the previous syntaxes. +## @var{outtype} can take the following values: ## -## If @var{x} is a matrix, compute the median value for each column and -## return them in a row vector. +## @table @asis +## @item "default" +## Output is of type double, unless the input is single in which case the output +## is of type single. +## +## @item "double" +## Output is of type double. ## -## If the optional @var{dim} argument is given, operate along this dimension. -## @seealso{mean, mode} +## @item "native". +## Output is of the same type as the input (@code{class (@var{x})}), unless the +## input is logical in which case the output is of type double. +## @end table +## +## The optional variable @var{nanflag} specifies whether to include or exclude +## NaN values from the calculation using any of the previously specified input +## argument combinations. The default value for @var{nanflag} is "includenan" +## which keeps NaN values in the calculation. To exclude NaN values set the +## value of @var{nanflag} to "omitnan". The output will still contain NaN +## values if @var{x} consists of all NaN values in the operating dimension. +## +## @seealso{mean, mode, movmedian} ## @end deftypefn -function y = median (x, dim) +function m = median (x, varargin) - if (nargin < 1) + if (nargin < 1 || nargin > 4) print_usage (); endif if (! (isnumeric (x) || islogical (x))) - error ("median: X must be a numeric vector or matrix"); + error ("median: X must be either numeric or logical."); + endif + + ## Set initial conditions + all_flag = 0; + omitnan = 0; + perm_flag = 0; + out_flag = 0; + vecdim_flag = 0; + dim = []; + + nvarg = numel (varargin); + varg_chars = cellfun ("ischar", varargin); + szx = sz_out = size (x); + ndx = ndims (x); + outtype = class (x); + + if (nvarg > 1 && ! varg_chars(2:end)) + ## Only first varargin can be numeric + print_usage (); + endif + + ## Process any other char arguments. + if (any (varg_chars)) + for idx = varargin(varg_chars) + switch (tolower (idx{:})) + case "all" + all_flag = 1; + + case "omitnan" + omitnan = 1; + + case "includenan" + omitnan = 0; + + case "native" + if (out_flag) + error ("median: only one OUTTYPE can be specified.") + endif + if (strcmp (outtype, "logical")) + outtype = "double"; + endif + out_flag = 1; + + case "default" + if (out_flag) + error ("median: only one OUTTYPE can be specified.") + endif + if (! strcmp (outtype, "single")) + outtype = "double"; + endif + out_flag = 1; + + case "double" + if (out_flag) + error ("median: only one OUTTYPE can be specified.") + endif + outtype = "double"; + out_flag = 1; + + otherwise + print_usage (); + endswitch + endfor + + varargin(varg_chars) = []; + nvarg = numel (varargin); + endif + + if (((nvarg == 1) && ! (isnumeric (varargin{1}))) || (nvarg > 1)) + ## After trimming char inputs should only be one numeric varargin left + print_usage (); + endif + + ## Process special cases for in/out size + if (nvarg > 0) + ## dim or vecdim provided + if (all_flag) + error ("median: 'all' cannot be used with DIM or VECDIM options."); + endif + + dim = varargin{1}; + vecdim_flag = ! isscalar (dim); + + if (! (isvector (dim) && (dim > 0)) || any (rem (dim, 1))) + error ("median: DIM must be a positive integer scalar or vector."); + endif + + ## Adjust sz_out, account for possible dim > ndx by appending singletons + sz_out(ndx + 1 : max (dim)) = 1; + sz_out(dim(dim <= ndx)) = 1; + szx(ndx + 1 : max (dim)) = 1; + + if (vecdim_flag) + ## vecdim - try to simplify first + dim = sort (dim); + if (! all (diff (dim))) + error ("median: VECDIM must contain non-repeating positive integers."); + endif + + ## dims > ndims(x) and dims only one element long don't affect median + sing_dim_x = find (szx != 1); + dim(dim > ndx | szx(dim) == 1) = []; + + if (isempty (dim)) + ## No dims left to process, return input as output + if (! strcmp (class (x), outtype)) + m = outtype_convert (x, outtype); + else + m = x; + endif + return; + elseif ((numel(dim) == numel(sing_dim_x)) ... + && unique ([dim, sing_dim_x]) == dim) + ## If DIMs cover all nonsingleton ndims(x) it's equivalent to "all" + ## (check lengths first to reduce unique overhead if not covered) + all_flag = true; + endif + endif + + else + ## Dim not provided. Determine scalar dimension + if (all_flag) + ## Special case 'all': Recast input as dim1 vector, process as normal + x = x(:); + szx = [numel(x), 1]; + dim = 1; + sz_out = [1 1]; + + elseif (isrow (x)) + ## Special case row vector: Avoid setting dim to 1. + dim = 2; + sz_out = [1, 1]; + + elseif (ndx == 2 && szx == [0, 0]) + ## Special case []: Do not apply sz_out(dim)=1 change + dim = 1; + sz_out = [1, 1]; + + else + ## General case: Set dim to first non-singleton, contract sz_out along dim + (dim = find (szx != 1, 1)) || (dim = 1); + sz_out(dim) = 1; + endif endif if (isempty (x)) - error ("median: X cannot be an empty matrix"); + ## Empty input - output NaN or class equivalent in pre-determined size + switch (outtype) + case {"double", "single"} + m = NaN (sz_out, outtype); + case ("logical") + m = false (sz_out); + otherwise + m = cast (NaN (sz_out), outtype); + endswitch + return; + endif + + if (szx(dim) == 1) + ## Operation along singleton dimension - nothing to do + if (! strcmp (class (x), outtype)) + m = outtype_convert (x, outtype); + else + m = x; + endif + return; + endif + + ## Permute dim to simplify all operations along dim1. At func. end ipermute. + + ## FIXME: for very large data sets, flattening all vecdim dimensions into dim1 + ## could hit index type limits + + if ((numel (dim) > 1) || (dim != 1 && ! isvector (x))) + perm_vect = 1 : ndx; + + if (! vecdim_flag) + ## Move dim to dim 1 + perm_vect([1, dim]) = [dim, 1]; + x = permute (x, perm_vect); + szx([1, dim]) = szx([dim, 1]); + dim = 1; + + else + ## Move vecdims to front + perm_vect(dim) = []; + perm_vect = [dim, perm_vect]; + x = permute (x, perm_vect); + + ## Reshape all vecdims into dim1 + num_dim = prod (szx(dim)); + szx(dim) = []; + szx = [ones(1, numel(dim)), szx]; + szx(1) = num_dim; + x = reshape (x, szx); + dim = 1; + endif + + perm_flag = true; + endif + + ## Find column locations of NaNs + nanfree = ! any (isnan (x), dim); + if (omitnan && nanfree(:)) + ## Don't use omitnan path if no NaNs are present. Prevents any data types + ## without a defined NaN from following slower omitnan codepath. + omitnan = 0; endif - nd = ndims (x); - sz = size (x); - if (nargin < 2) - ## Find the first non-singleton dimension. - (dim = find (sz > 1, 1)) || (dim = 1); + x = sort (x, dim); # Note: pushes any NaN's to end for omitnan compatability + + if (omitnan) + ## Ignore any NaN's in data. Each operating vector might have a + ## different number of non-NaN data points. + + if (isvector (x)) + ## Checks above ensure either dim1 or dim2 vector + x = x(! isnan (x)); + n = numel (x); + k = floor ((n + 1) / 2); + if (mod (n, 2)) + ## odd + m = x(k); + else + ## even + m = (x(k) + x(k + 1)) / 2; + endif + + else + n = sum (! isnan (x), 1); + k = floor ((n + 1) ./ 2); + m_idx_odd = mod (n, 2) & n; + m_idx_even = (! m_idx_odd) & n; + + m = NaN ([1, szx(2 : end)]); + + if (ndims (x) > 2) + szx = [szx(1), prod(szx(2 : end))]; + endif + + ## Grab kth value, k possibly different for each column + if (any (m_idx_odd(:))) + x_idx_odd = sub2ind (szx, k(m_idx_odd), find (m_idx_odd)); + m(m_idx_odd) = x(x_idx_odd); + endif + if (any (m_idx_even(:))) + k_even = k(m_idx_even)(:); + x_idx_even = sub2ind (szx, [k_even, k_even + 1], ... + (find (m_idx_even))(:, [1, 1])); + m(m_idx_even) = sum (x(x_idx_even), 2) / 2; + endif + endif + else - if (! (isscalar (dim) && dim == fix (dim) && dim > 0)) - error ("median: DIM must be an integer and a valid dimension"); + ## No "omitnan". All 'vectors' uniform length. + ## All types without a NaN value will use this path. + if (all (! nanfree)) + m = NaN (sz_out); + + else + if (isvector (x)) + n = numel (x); + k = floor ((n + 1) / 2); + + m = x(k); + if (! mod (n, 2)) + ## Even + if (any (isa (x, "integer"))) + ## avoid int overflow issues + m2 = x(k + 1); + if (sign (m) != sign (m2)) + m += m2; + m /= 2; + else + m += (m2 - m) / 2; + endif + else + m += (x(k + 1) - m) / 2; + endif + endif + + else + ## Nonvector, all operations were permuted to be along dim 1 + n = szx(1); + k = floor ((n + 1) / 2); + + if (isfloat (x)) + m = NaN ([1, szx(2 : end)]); + else + m = zeros ([1, szx(2 : end)], outtype); + endif + + if (! mod (n, 2)) + ## Even + if (any (isa (x, "integer"))) + ## avoid int overflow issues + + ## Use flattened index to simplify n-D operations + m(1, :) = x(k, :); + m2 = x(k + 1, :); + + samesign = prod (sign ([m(1, :); m2]), 1) == 1; + m(1, :) = samesign .* m(1, :) + ... + (m2 + !samesign .* m(1, :) - samesign .* m(1, :)) / 2; + + else + m(nanfree) = (x(k, nanfree) + x(k + 1, nanfree)) / 2; + endif + else + ## Odd. Use flattened index to simplify n-D operations + m(nanfree) = x(k, nanfree); + endif + endif endif endif - n = size (x, dim); - k = floor ((n+1) / 2); - if (mod (n, 2) == 1) - y = nth_element (x, k, dim); - else - y = sum (nth_element (x, k:k+1, dim), dim, "native") / 2; - if (islogical (x)) - y = logical (y); - endif + if (perm_flag) + ## Inverse permute back to correct dimensions + m = ipermute (m, perm_vect); endif - ## Inject NaNs where needed, to be consistent with Matlab. - if (isfloat (x)) - y(any (isnan (x), dim)) = NaN; + + ## Convert output type as requested + if (! strcmp (class (m), outtype)) + m = outtype_convert (m, outtype); endif endfunction +function m = outtype_convert (m, outtype) + switch (outtype) + case "single" + m = single (m); + case "double" + m = double (m); + otherwise + m = cast (m, outtype); + endswitch +endfunction + +%!assert (median (1), 1) +%!assert (median ([1,2,3]), 2) +%!assert (median ([1,2,3]'), 2) +%!assert (median (cat(3,3,1,2)), 2) +%!assert (median ([3,1,2]), 2) +%!assert (median ([2,4,6,8]), 5) +%!assert (median ([8,2,6,4]), 5) +%!assert (median (single ([1,2,3])), single (2)) +%!assert (median ([1,2], 3), [1,2]) %!test %! x = [1, 2, 3, 4, 5, 6]; @@ -111,12 +464,159 @@ %! %! assert (median (x) == median (x2) && median (x) == 3.5); %! assert (median (y) == median (y2) && median (y) == 4); -%! assert (median ([x2, 2*x2]), [3.5, 7]); -%! assert (median ([y2, 3*y2]), [4, 12]); +%! assert (median ([x2, 2 * x2]), [3.5, 7]); +%! assert (median ([y2, 3 * y2]), [4, 12]); + +## Test outtype option +%!test +%! in = [1 2 3]; +%! out = 2; +%! assert (median (in, "default"), median (in)); +%! assert (median (in, "default"), out); +%!test +%! in = single ([1 2 3]); +%! out = 2; +%! assert (median (in, "default"), single (median (in))); +%! assert (median (in, "default"), single (out)); +%! assert (median (in, "double"), double (out)); +%! assert (median (in, "native"), single (out)); +%!test +%! in = uint8 ([1 2 3]); +%! out = 2; +%! assert (median (in, "default"), double (median (in))); +%! assert (median (in, "default"), double (out)); +%! assert (median (in, "double"), out); +%! assert (median (in, "native"), uint8 (out)); +%!test +%! in = logical ([1 0 1]); +%! out = 1; +%! assert (median (in, "default"), double (median (in))); +%! assert (median (in, "default"), double (out)); +%! assert (median (in, "double"), double (out)); +%! assert (median (in, "native"), double (out)); + +## Test single input and optional arguments "all", DIM, "omitnan") +%!test +%! x = repmat ([2 2.1 2.2 2 NaN; 3 1 2 NaN 5; 1 1.1 1.4 5 3], [1, 1, 4]); +%! y = repmat ([2 1.1 2 NaN NaN], [1, 1, 4]); +%! assert (median (x), y); +%! assert (median (x, 1), y); +%! y = repmat ([2 1.1 2 3.5 4], [1, 1, 4]); +%! assert (median (x, "omitnan"), y); +%! assert (median (x, 1, "omitnan"), y); +%! y = repmat ([2.05; 2.5; 1.4], [1, 1, 4]); +%! assert (median (x, 2, "omitnan"), y); +%! y = repmat ([NaN; NaN; 1.4], [1, 1, 4]); +%! assert (median (x, 2), y); +%! assert (median (x, "all"), NaN); +%! assert (median (x, "all", "omitnan"), 2); +%!assert (median (cat (3, 3, 1, NaN, 2), "omitnan"), 2) +%!assert (median (cat (3, 3, 1, NaN, 2), 3, "omitnan"), 2) + +# Test boolean input +%!test +%! assert (median (true, "all"), logical (1)); +%! assert (median (false), logical (0)); +%! assert (median ([true false true]), true); +%! assert (median ([true false true], 2), true); +%! assert (median ([true false true], 1), logical ([1 0 1])); +%! assert (median ([true false NaN], 1), [1 0 NaN]); +%! assert (median ([true false NaN], 2), NaN); +%! assert (median ([true false NaN], 2, "omitnan"), 0.5); +%! assert (median ([true false NaN], 2, "omitnan", "native"), double(0.5)); + +## Test dimension indexing with vecdim in n-dimensional arrays +%!test +%! x = repmat ([1:20;6:25], [5 2 6 3]); +%! assert (size (median (x, [3 2])), [10 1 1 3]); +%! assert (size (median (x, [1 2])), [1 1 6 3]); +%! assert (size (median (x, [1 2 4])), [1 1 6]); +%! assert (size (median (x, [1 4 3])), [1 40]); +%! assert (size (median (x, [1 2 3 4])), [1 1]); + +## Test exceeding dimensions +%!assert (median (ones (2,2), 3), ones (2,2)); +%!assert (median (ones (2,2,2), 99), ones (2,2,2)); +%!assert (median (magic (3), 3), magic (3)); +%!assert (median (magic (3), [1 3]), [4, 5, 6]); +%!assert (median (magic (3), [1 99]), [4, 5, 6]); -%!assert (median (single ([1,2,3])), single (2)) +## Test results with vecdim in n-dimensional arrays and "omitnan" +%!test +%! x = repmat ([2 2.1 2.2 2 NaN; 3 1 2 NaN 5; 1 1.1 1.4 5 3], [1, 1, 4]); +%! assert (median (x, [3 2]), [NaN NaN 1.4]'); +%! assert (median (x, [3 2], "omitnan"), [2.05 2.5 1.4]'); +%! assert (median (x, [1 3]), [2 1.1 2 NaN NaN]); +%! assert (median (x, [1 3], "omitnan"), [2 1.1 2 3.5 4]); + +## Test empty, NaN, Inf inputs +%!assert (median (NaN), NaN) +%!assert (median (NaN, "omitnan"), NaN) +%!assert (median (NaN (2)), [NaN NaN]) +%!assert (median (NaN (2), "omitnan"), [NaN NaN]) +%!assert (median ([1 NaN 3]), NaN) +%!assert (median ([1 NaN 3], 1), [1 NaN 3]) +%!assert (median ([1 NaN 3], 2), NaN) +%!assert (median ([1 NaN 3]'), NaN) +%!assert (median ([1 NaN 3]', 1), NaN) +%!assert (median ([1 NaN 3]', 2), [1; NaN; 3]) +%!assert (median ([1 NaN 3], "omitnan"), 2) +%!assert (median ([1 NaN 3]', "omitnan"), 2) +%!assert (median ([1 NaN 3], 1, "omitnan"), [1 NaN 3]) +%!assert (median ([1 NaN 3], 2, "omitnan"), 2) +%!assert (median ([1 NaN 3]', 1, "omitnan"), 2) +%!assert (median ([1 NaN 3]', 2, "omitnan"), [1; NaN; 3]) +%!assert (median ([1 2 NaN 3]), NaN) +%!assert (median ([1 2 NaN 3], "omitnan"), 2) %!assert (median ([1,2,NaN;4,5,6;NaN,8,9]), [NaN, 5, NaN]) -%!assert (median ([1,2], 3), [1,2]) +%!assert (median ([1 2 ; NaN 4]), [NaN 3]) +%!assert (median ([1 2 ; NaN 4], "omitnan"), [1 3]) +%!assert (median ([1 2 ; NaN 4], 1, "omitnan"), [1 3]) +%!assert (median ([1 2 ; NaN 4], 2, "omitnan"), [1.5; 4], eps) +%!assert (median ([1 2 ; NaN 4], 3, "omitnan"), [1 2 ; NaN 4]) +%!assert (median ([NaN 2 ; NaN 4]), [NaN 3]) +%!assert (median ([NaN 2 ; NaN 4], "omitnan"), [NaN 3]) +%!assert (median (ones (1, 0, 3)), NaN (1, 1, 3)) + +%!assert (median (NaN("single")), NaN("single")); +%!assert (median (NaN("single"), "omitnan"), NaN("single")); +%!assert (median (NaN("single"), "double"), NaN("double")); +%!assert (median (single([1 2 ; NaN 4])), single([NaN 3])); +%!assert (median (single([1 2 ; NaN 4]), "double"), double([NaN 3])); +%!assert (median (single([1 2 ; NaN 4]), "omitnan"), single([1 3])); +%!assert (median (single([1 2 ; NaN 4]), "omitnan", "double"), double([1 3])); +%!assert (median (single([NaN 2 ; NaN 4]), "double"), double([NaN 3])); +%!assert (median (single([NaN 2 ; NaN 4]), "omitnan"), single([NaN 3])); +%!assert (median (single([NaN 2 ; NaN 4]), "omitnan", "double"), double([NaN 3])); + +%!assert (median (Inf), Inf); +%!assert (median (-Inf), -Inf); +%!assert (median ([-Inf Inf]), NaN); +%!assert (median ([3 Inf]), Inf); +%!assert (median ([3 4 Inf]), 4); +%!assert (median ([Inf 3 4]), 4); +%!assert (median ([Inf 3 Inf]), Inf); + +%!assert (median ([]), NaN); +%!assert (median (ones(1,0)), NaN); +%!assert (median (ones(0,1)), NaN); +%!assert (median ([], 1), NaN(1,0)); +%!assert (median ([], 2), NaN(0,1)); +%!assert (median ([], 3), NaN(0,0)); +%!assert (median (ones(1,0), 1), NaN(1,0)); +%!assert (median (ones(1,0), 2), NaN(1,1)); +%!assert (median (ones(1,0), 3), NaN(1,0)); +%!assert (median (ones(0,1), 1), NaN(1,1)); +%!assert (median (ones(0,1), 2), NaN(0,1)); +%!assert (median (ones(0,1), 3), NaN(0,1)); +%!assert (median (ones(0,1,0,1), 1), NaN(1,1,0)); +%!assert (median (ones(0,1,0,1), 2), NaN(0,1,0)); +%!assert (median (ones(0,1,0,1), 3), NaN(0,1,1)); +%!assert (median (ones(0,1,0,1), 4), NaN(0,1,0)); + +## Test complex inputs (should sort by abs(a)) +%!assert (median([1 3 3i 2 1i]), 2) +%!assert (median([1 2 4i; 3 2i 4]), [2, 1+1i, 2+2i]) ## Test multidimensional arrays %!shared a, b, x, y @@ -129,18 +629,80 @@ %! y = sort (b, 3); %!assert <*35679> (median (a, 4), x(:, :, :, 3)) %!assert <*35679> (median (b, 3), (y(:, :, 3, :) + y(:, :, 4, :))/2) +%!shared ## Clear shared to prevent variable echo for any later test failures ## Test non-floating point types %!assert (median ([true, false]), true) +%!assert (median (logical ([])), false) %!assert (median (uint8 ([1, 3])), uint8 (2)) +%!assert (median (uint8 ([])), uint8 (NaN)) +%!assert (median (uint8 ([NaN 10])), uint8 (5)) %!assert (median (int8 ([1, 3, 4])), int8 (3)) +%!assert (median (int8 ([])), int8 (NaN)) %!assert (median (single ([1, 3, 4])), single (3)) %!assert (median (single ([1, 3, NaN])), single (NaN)) +## Test same sign int overflow when getting mean of even number of values +%!assert <*54567> (median (uint8 ([253, 255])), uint8 (254)) +%!assert <*54567> (median (uint8 ([253, 254])), uint8 (254)) +%!assert <*54567> (median (int8 ([127, 126, 125, 124; 1 3 5 9])), ... +%! int8 ([64 65 65 67])) +%!assert <*54567> (median (int8 ([127, 126, 125, 124; 1 3 5 9]), 2), ... +%! int8 ([126; 4])) +%!assert <*54567> (median (int64 ([intmax("int64"), intmax("int64")-2])), ... +%! intmax ("int64") - 1) +%!assert <*54567> (median ( ... +%! int64 ([intmax("int64"), intmax("int64")-2; 1 2]), 2), ... +%! int64([intmax("int64") - 1; 2])) +%!assert <*54567> (median (uint64 ([intmax("uint64"), intmax("uint64")-2])), ... +%! intmax ("uint64") - 1) +%!assert <*54567> (median ( ... +%! uint64 ([intmax("uint64"), intmax("uint64")-2; 1 2]), 2), ... +%! uint64([intmax("uint64") - 1; 2])) + +## Test opposite sign int overflow when getting mean of even number of values +%!assert <*54567> (median (... +%! [intmin('int8') intmin('int8')+5 intmax('int8')-5 intmax('int8')]), ... +%! int8(-1)) +%!assert <*54567> (median ([int8([1 2 3 4]); ... +%! intmin('int8') intmin('int8')+5 intmax('int8')-5 intmax('int8')], 2), ... +%! int8([3;-1])) +%!assert <*54567> (median (... +%! [intmin('int64') intmin('int64')+5 intmax('int64')-5 intmax('int64')]), ... +%! int64(-1)) +%!assert <*54567> (median ([int64([1 2 3 4]); ... +%! intmin('int64') intmin('int64')+5 intmax('int64')-5 intmax('int64')], 2), ... +%! int64([3;-1])) + +## Test int accuracy loss doing mean of close int64/uint64 values as double +%!assert <*54567> (median ([intmax("uint64"), intmax("uint64")-2]), ... +%! intmax("uint64")-1) +%!assert <*54567> (median ([intmax("uint64"), intmax("uint64")-2], "default"), ... +%! double(intmax("uint64")-1)) +%!assert <*54567> (median ([intmax("uint64"), intmax("uint64")-2], "double"), ... +%! double(intmax("uint64")-1)) +%!assert <*54567> (median ([intmax("uint64"), intmax("uint64")-2], "native"), ... +%! intmax("uint64")-1) + +## Test input case insensitivity +%!assert (median ([1 2 3], "aLL"), 2); +%!assert (median ([1 2 3], "OmitNan"), 2); +%!assert (median ([1 2 3], "DOUBle"), 2); + ## Test input validation %!error <Invalid call> median () -%!error <X must be a numeric> median ({1:5}) -%!error <X cannot be an empty matrix> median ([]) -%!error <DIM must be an integer> median (1, ones (2,2)) -%!error <DIM must be an integer> median (1, 1.5) -%!error <DIM must be .* a valid dimension> median (1, 0) +%!error <Invalid call> median (1, 2, 3) +%!error <Invalid call> median (1, 2, 3, 4) +%!error <Invalid call> median (1, "all", 3) +%!error <Invalid call> median (1, "b") +%!error <Invalid call> median (1, 1, "foo") +%!error <'all' cannot be used with> median (1, 3, "all") +%!error <'all' cannot be used with> median (1, [2 3], "all") +%!error <X must be either numeric or logical> median ({1:5}) +%!error <X must be either numeric or logical> median ("char") +%!error <only one OUTTYPE can be specified> median(1, "double", "native") +%!error <DIM must be a positive integer> median (1, ones (2,2)) +%!error <DIM must be a positive integer> median (1, 1.5) +%!error <DIM must be a positive integer> median (1, 0) +%!error <DIM must be a positive integer> median ([1 2 3], [-1 1]) +%!error <VECDIM must contain non-repeating> median(1, [1 2 2])
--- a/scripts/statistics/std.m Tue Feb 28 15:45:03 2023 -0500 +++ b/scripts/statistics/std.m Wed Mar 01 00:38:54 2023 -0500 @@ -24,84 +24,104 @@ ######################################################################## ## -*- texinfo -*- -## @deftypefn {} {@var{y} =} std (@var{x}) -## @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{}) +## @deftypefn {} {@var{s} =} std (@var{x}) +## @deftypefnx {} {@var{s} =} std (@var{x}, @var{w}) +## @deftypefnx {} {@var{s} =} std (@var{x}, @var{w}, @var{dim}) +## @deftypefnx {} {@var{s} =} std (@var{x}, @var{w}, @var{vecdim}) +## @deftypefnx {} {@var{s} =} std (@var{x}, @var{w}, @qcode{"ALL"}) +## @deftypefnx {} {@var{s} =} std (@dots{}, @var{nanflag}) +## @deftypefnx {} {[@var{s}, @var{m}] =} std (@dots{}) ## Compute the standard deviation of the elements of the vector @var{x}. ## ## The standard deviation is defined as ## @tex -## $$ -## {\rm std} (x) = \sigma = \sqrt{{\sum_{i=1}^N (x_i - \bar{x})^2 \over N - 1}} -## $$ -## where $\bar{x}$ is the mean value of @var{x} and $N$ is the number of elements of @var{x}. +## $$ {\rm std}(x) = \sqrt{{1\over N-1} \sum_{i=1}^N (x_i - \bar x )^2} $$ +## where $\bar{x}$ is the mean value of @var{x} and $N$ is the number of +## elements of @var{x}. ## @end tex ## @ifnottex ## ## @example ## @group -## std (@var{x}) = sqrt ( 1/(N-1) SUM_i (@var{x}(i) - mean(@var{x}))^2 ) +## std (@var{x}) = sqrt ((1 / (N-1)) * SUM_i ((@var{x}(i) - mean(@var{x}))^2)) ## @end group ## @end example ## ## @noindent -## where @math{N} is the number of elements of the @var{x} vector. +## where @math{N} is the number of elements of @var{x}. ## @end ifnottex ## -## If @var{x} is an array, compute the standard deviation for each column and -## return them in a row vector (or for an n-D array, the result is returned as -## an array of dimension 1 x n x m x @dots{}). +## If @var{x} is an array, compute the standard deviation along the first +## non-singleton dimensions of @var{x}. ## ## The optional argument @var{w} determines the weighting scheme to use. Valid ## values are: ## ## @table @asis ## @item 0 [default]: -## Normalize with @math{N-1}. This provides the square root of the best -## unbiased estimator of the variance. +## Normalize with @math{N-1} (population standard deviation). This provides the +## square root of the best unbiased estimator of the standard deviation. ## ## @item 1: -## Normalize with @math{N}@. This provides the square root of the second -## moment around the mean. +## Normalize with @math{N} (sample standard deviation). This provides the +## square root of the second moment around the mean. ## ## @item a vector: -## Compute the weighted standard deviation with non-negative scalar weights. -## The length of @var{w} must equal the size of @var{x} along dimension -## @var{dim}. +## Compute the weighted standard deviation with non-negative weights. +## The length of @var{w} must equal the size of @var{x} in the operating +## dimension. NaN values are permitted in @var{w}, will be multiplied with the +## associated values in @var{x}, and can be excluded by the @var{nanflag} +## option. +## +## @item an array: +## Similar to vector weights, but @var{w} must be the same size as @var{x}. If +## the operating dimension is supplied as @var{vecdim} or "all" and @var{w} is +## not a scalar, @var{w} must be an same-sized array. ## @end table ## -## If @math{N} is equal to 1 the value of @var{W} is ignored and normalization -## by @math{N} is used. +## Note: @var{w} must always be specified before specifying any of the following +## dimension options. To use the default value for @var{w} you may pass an empty +## input argument []. ## ## The optional variable @var{dim} forces @code{std} to operate over the -## specified dimension(s). @var{dim} can either be a scalar dimension or a -## vector of non-repeating dimensions. Dimensions must be positive integers, -## and the standard deviation is calculated over the array slice defined by -## @var{dim}. +## specified dimension, which must be a positive integer-valued number. +## Specifying any singleton dimension in @var{x}, including any dimension +## exceeding @code{ndims (@var{x})}, will result in a standard deviation of 0. +## +## Specifying the dimensions as @var{vecdim}, a vector of non-repeating +## dimensions, will return the standard deviation calculated over the array +## slice defined by @var{vecdim}. If @var{vecdim} indexes all dimensions of +## @var{x}, then it is equivalent to the option @qcode{"all"}. Any dimension in +## @var{vecdim} greater than @code{ndims (@var{x})} is ignored. ## -## Specifying dimension @qcode{"all"} will force @code{std} to operate on all -## elements of @var{x}, and is equivalent to @code{std (@var{x}(:))}. +## Specifying the dimension as @qcode{"all"} will force @code{std} to operate on +## all 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. +## The optional variable @var{nanflag} specifies whether to include or exclude +## NaN values from the calculation using any of the previously specified input +## argument combinations. The default value for @var{nanflag} is "includenan" +## which keeps NaN values in the calculation. To exclude NaN values set the +## value of @var{nanflag} to "omitnan". The output will still contain NaN +## values if @var{x} consists of all NaN values in the operating dimension. ## -## The optional second output variable @var{mu} contains the mean or weighted -## mean used to calculate @var{y}, and will be the same size as @var{y}. +## The optional second output variable @var{mu} contains the mean of the +## elements of @var{x} used to calculate the standard deviation. If @var{v} is +## the weighted standard deviation, then @var{m} is also the weighted mean. +## ## @seealso{var, bounds, mad, range, iqr, mean, median} ## @end deftypefn -function [y, mu] = std (varargin) +function [s, m] = std (varargin) if (nargin < 1) print_usage (); endif if (nargout < 2) - y = sqrt (var (varargin{:})); + s = sqrt (var (varargin{:})); else - [y, mu] = var (varargin{:}); - y = sqrt (y); + [s, m] = var (varargin{:}); + s = sqrt (s); endif endfunction
--- a/scripts/statistics/var.m Tue Feb 28 15:45:03 2023 -0500 +++ b/scripts/statistics/var.m Wed Mar 01 00:38:54 2023 -0500 @@ -1,6 +1,6 @@ ######################################################################## ## -## Copyright (C) 1995-2023 The Octave Project Developers +## Copyright (C) 1996-2023 The Octave Project Developers ## ## See the file COPYRIGHT.md in the top-level directory of this ## distribution or <https://octave.org/copyright/>. @@ -27,214 +27,542 @@ ## @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}, @var{vecdim}) ## @deftypefnx {} {@var{v} =} var (@var{x}, @var{w}, @qcode{"ALL"}) +## @deftypefnx {} {@var{v} =} var (@dots{}, @var{nanflag}) ## @deftypefnx {} {[@var{v}, @var{m}] =} var (@dots{}) ## Compute the variance of the elements of the vector @var{x}. ## ## The variance is defined as ## @tex -## $$ -## {\rm var} (x) = \sigma^2 = {\sum_{i=1}^N (x_i - \bar{x})^2 \over N - 1} -## $$ +## $$ {\rm var}(x) = {1\over N-1} \sum_{i=1}^N (x_i - \bar x )^2 $$ ## where $\bar{x}$ is the mean value of @var{x} and $N$ is the number of ## elements of @var{x}. -## ## @end tex ## @ifnottex ## ## @example ## @group -## var (@var{x}) = 1/(N-1) SUM_i (@var{x}(i) - mean(@var{x}))^2 +## var (@var{x}) = (1 / (N-1)) * SUM_i ((@var{x}(i) - mean(@var{x}))^2) ## @end group ## @end example ## ## @noindent -## where @math{N} is the length of the @var{x} vector. -## +## where @math{N} is the number of elements of @var{x}. ## @end ifnottex -## If @var{x} is an array, compute the variance for each column and return them -## in a row vector (or for an n-D array, the result is returned as an array of -## dimension 1 x n x m x @dots{}). +## +## If @var{x} is an array, compute the variance along the first non-singleton +## dimensions of @var{x}. ## ## The optional argument @var{w} determines the weighting scheme to use. Valid -## values are +## values are: ## ## @table @asis ## @item 0 [default]: -## Normalize with @math{N-1}. This provides the square root of the best -## unbiased estimator of the variance. +## Normalize with @math{N-1} (population variance). This provides the square +## root of the best unbiased estimator of the variance. ## ## @item 1: -## Normalize with @math{N}@. This provides the square root of the second -## moment around the mean. +## Normalize with @math{N} (sample variance). This provides the square root of +## the second moment around the mean. ## ## @item a vector: -## Compute the weighted variance with non-negative scalar weights. The length -## of @var{w} must equal the size of @var{x} along dimension @var{dim}. +## Compute the weighted variance with non-negative weights. The length of +## @var{w} must equal the size of @var{x} in the operating dimension. NaN +## values are permitted in @var{w}, will be multiplied with the associated +## values in @var{x}, and can be excluded by the @var{nanflag} option. +## +## @item an array: +## Similar to vector weights, but @var{w} must be the same size as @var{x}. If +## the operating dimension is supplied as @var{vecdim} or "all" and @var{w} is +## not a scalar, @var{w} must be an same-sized array. ## @end table ## -## If @math{N} is equal to 1 the value of @var{W} is ignored and normalization -## by @math{N} is used. +## Note: @var{w} must always be specified before specifying any of the following +## dimension options. To use the default value for @var{w} you may pass an empty +## input argument []. ## ## The optional variable @var{dim} forces @code{var} to operate over the -## specified dimension(s). @var{dim} can either be a scalar dimension or a -## vector of non-repeating dimensions. Dimensions must be positive integers, -## and the variance is calculated over the array slice defined by @var{dim}. +## specified dimension, which must be a positive integer-valued number. +## Specifying any singleton dimension in @var{x}, including any dimension +## exceeding @code{ndims (@var{x})}, will result in a variance of 0. +## +## Specifying the dimensions as @var{vecdim}, a vector of non-repeating +## dimensions, will return the variance calculated over the array slice defined +## by @var{vecdim}. If @var{vecdim} indexes all dimensions of @var{x}, then it +## is equivalent to the option @qcode{"all"}. Any dimension in @var{vecdim} +## greater than @code{ndims (@var{x})} is ignored. +## +## Specifying the dimension as @qcode{"all"} will force @code{var} to operate on +## all elements of @var{x}, and is equivalent to @code{var (@var{x}(:))}. ## -## Specifying dimension @qcode{"all"} will force @code{var} to operate on all -## elements of @var{x}, and is equivalent to @code{var (@var{x}(:))}. +## The optional variable @var{nanflag} specifies whether to include or exclude +## NaN values from the calculation using any of the previously specified input +## argument combinations. The default value for @var{nanflag} is "includenan" +## which keeps NaN values in the calculation. To exclude NaN values set the +## value of @var{nanflag} to "omitnan". The output will still contain NaN +## values if @var{x} consists of all NaN values in the operating dimension. ## -## When @var{dim} is a vector or @qcode{"all"}, @var{w} must be either 0 or 1. +## The optional second output variable @var{mu} contains the mean of the +## elements of @var{x} used to calculate the variance. If @var{v} is the +## weighted variance, then @var{m} is also the weighted mean. ## -## The optional second output variable @var{mu} contains the mean or weighted -## mean used to calculate @var{v}, and will be the same size as @var{v}. -## @seealso{cov, std, skewness, kurtosis, moment} +## @seealso{std, mean, cov, skewness, kurtosis, moment} ## @end deftypefn -function [v, mu] = var (x, w = 0, dim = []) +function [v, m] = var (x, varargin) + + if (nargin < 1 || nargin > 4) + print_usage (); + endif - if (nargin < 1) + ## initialize variables + all_flag = false; + omitnan = false; + nvarg = numel (varargin); + varg_chars = cellfun ('ischar', varargin); + + ## Check all char arguments. + if (nvarg == 3 && ! varg_chars(3)) print_usage (); endif - if (! (isnumeric (x) || islogical (x))) - error ("var: X must be a numeric vector or matrix"); + if (any (varg_chars)) + for i = varargin(varg_chars) + switch (lower (i{:})) + case "all" + all_flag = true; + case "omitnan" + omitnan = true; + case "includenan" + omitnan = false; + otherwise + print_usage (); + endswitch + endfor + varargin(varg_chars) = []; + nvarg = numel (varargin); endif - if (isempty (w)) - w = 0; + # FIXME: when sparse can use broadcast ops, remove sparse checks and hacks + sprs_x = issparse (x); + w = 0; + weighted = false; # true if weight vector/array used + vecdim = []; + vecempty = true; + vecdim_scalar_vector = [false, false]; # [false, false] for empty vecdim + szx = size (x); + ndx = ndims (x); + + ## Check numeric arguments + if (! (isnumeric (x))) + error ("var: X must be a numeric vector or matrix."); + endif + if (isa (x, "single")) + outtype = "single"; + else + outtype = "double"; endif - nd = ndims (x); - sz = size (x); - emptydimflag = false; - - if (isempty (dim)) - emptydimflag = true; # Compatibility hack for empty x, ndims==2 - - ## Find the first non-singleton dimension. - (dim = find (sz != 1, 1)) || (dim = 1); + if (nvarg > 0) + if (nvarg > 2 || any (! cellfun ('isnumeric', varargin))) + print_usage (); + endif + ## Process weight input + if (any (varargin{1} < 0)) + error ("var: weights must not contain any negative values."); + endif + if (isscalar (varargin{1})) + w = varargin{1}; + if (! (w == 0 || w == 1) && ! isscalar (x)) + error ("var: normalization scalar must be either 0 or 1."); + endif + elseif (numel (varargin{1}) > 1) + weights = varargin{1}; + weighted = true; + endif + if (nvarg > 1) + ## Process dimension input + vecdim = varargin{2}; + if (! (vecempty = isempty (vecdim))) + ## Check for empty vecdim, won't change vsv if nonzero size empty + vecdim_scalar_vector = [isscalar(vecdim), isvector(vecdim)]; + endif + if (! (vecdim_scalar_vector(2) && all (vecdim > 0)) ... + || any (rem (vecdim, 1))) + error ("var: DIM must be a positive integer scalar or vector."); + endif + if (vecdim_scalar_vector(1) && vecdim > ndx && ! isempty (x)) + ## Scalar dimension larger than ndims(x), variance of any single number + ## is zero, except for inf, NaN, and empty values of x. + v = zeros (szx, outtype); + vn = ! isfinite (x); + v(vn) = NaN; + m = x; + return; + endif + if (vecdim_scalar_vector == [0 1] && (! all (diff (sort (vecdim))))) + error ("var: VECDIM must contain non-repeating positive integers."); + endif + endif + endif - else - if (isscalar (dim)) - if (dim < 1 || dim != fix (dim)) - error ("var: DIM must be a positive integer scalar, vector, or 'all'"); - endif - elseif (isnumeric (dim)) - if (! isvector (dim) && all (dim > 0) && all (rem (dim, 1) == 0)) - error ("var: DIM must be a positive integer scalar, vector, or 'all'"); - endif - if (dim != unique (dim, "stable")) - error ("var: vector DIM must contain non-repeating positive integers"); - endif - if (! isscalar (w)) - error ("var: W must be either 0 or 1 when DIM is a vector"); + ## Check for conflicting input arguments + if (all_flag && ! vecempty) + error ("var: 'all' flag cannot be used with DIM or VECDIM options."); + endif + if (weighted) + if (all_flag) + if (isvector (weights)) + if (numel (weights) != numel (x)) + error ("var: weight vector element count does not match X."); + endif + elseif (! (isequal (size (weights), szx))) + error ("var: weight matrix or array does not match X in size."); endif - ## Reshape X to compute the variance over an array slice - if (iscolumn (dim)) - dim = dim.'; + elseif (vecempty) + dim = find (szx > 1, 1); + if length (dim) == 0 + dim = 1; + endif + if (isvector (weights)) + if (numel (weights) != szx(dim)) + error (["var: weight vector length does not match operating ", ... + "dimension."]); + endif + elseif (! isequal (size (weights), szx)) + error ("var: weight matrix or array does not match X in size."); + endif + elseif (vecdim_scalar_vector(1)) + if (isvector (weights)) + if (numel (weights) != szx(vecdim)) + error (["var: weight vector length does not match operating ", ... + "dimension."]); + endif + elseif (! isequal (size (weights), szx)) + error ("var: weight matrix or array does not match X in size."); endif - collapsed_dims = dim; - dim = dim(end); - - ## Permute X to cluster the dimensions to collapse - max_dim = max ([nd, collapsed_dims]); - perm_start = perm_end = [1:max_dim]; - perm_start(dim:end) = []; - perm_start(ismember (perm_start, collapsed_dims)) = []; - perm_end(1:dim) = []; - perm_end(ismember (perm_end, collapsed_dims)) = []; - perm = [perm_start, collapsed_dims, perm_end]; - - x = permute (x, perm); + elseif (vecdim_scalar_vector(2) && ! (isequal (size (weights), szx))) + error ("var: weight matrix or array does not match X in size."); + endif + endif - ## Collapse the given dimensions - newshape = ones (1, max_dim); - newshape(1:nd) = sz; - newshape(collapsed_dims(1:(end-1))) = 1; - newshape(dim) = prod (sz(collapsed_dims)); - - ## New X with collapsed dimensions - x = reshape (x, newshape); - - elseif (ischar (dim) && strcmpi (dim, "all")) - if (! isscalar (w)) - error ("var: W must be either 0 or 1 when using 'all' as dimension"); + ## Force output for X being empty or scalar + if (isempty (x)) + if (vecempty && (ndx == 2 || all ((szx) == 0))) + v = NaN (outtype); + if (nargout > 1) + m = NaN (outtype); endif - - ## "all" equates to collapsing all elements to a single vector - x = x(:); - dim = 1; - sz = size (x); - else - error ("var: DIM must be a positive integer scalar, vector, or 'all'"); + return; + endif + if (vecdim_scalar_vector(1)) + szx(vecdim) = 1; + v = NaN (szx, outtype); + if (nargout > 1) + m = NaN (szx, outtype); + endif + return; endif endif - n = size (x, dim); - if (! isvector (w) || ! isnumeric (w) - || (isvector (w) && any (w < 0)) || - (isscalar (w) && ((w != 0 && w != 1) && (n != 1)))) - error ("var: W must be 0, 1, or a vector of positive integers"); + if (isscalar (x)) + if (isfinite (x)) + v = zeros (outtype); + else + v = NaN (outtype); + endif + if (nargout > 1) + m = x; + endif + return; endif - if (isempty (x)) - ## Empty matrix special case - if (emptydimflag && nd == 2 && all (sz == [0, 0])) - v = NaN; - mu = NaN; - else - sz(dim) = 1; - v = NaN (sz); - mu = NaN (sz); - endif - elseif (n == 1) - ## Scalar special case - if (! isscalar (w)) - error (["var: the length of W must be equal to the size of X " ... - "in the dimension along which variance is calculated"]); - endif - if (isa (x, "single")) - v = zeros (sz, "single"); - mu = x; - else - v = zeros (sz); - mu = x; - endif - v(isnan (x) | isinf (x)) = NaN; - else - ## Regular algorithm - if (isscalar (w)) - v = sumsq (center (x, dim), dim) / (n - 1 + w); - if (nargout == 2) - mu = mean (x, dim); + if (nvarg == 0) + ## Only numeric input argument, no dimensions or weights. + if (all_flag) + x = x(:); + if (omitnan) + x = x(! isnan (x)); + endif + n = length (x); + m = sum (x) ./ n; + v = sum (abs (x - m) .^ 2) ./ (n - 1 + w); + if (n == 1) + v = 0; endif else - ## Weighted variance - if (numel (w) != n) - error (["var: the length of W must be equal to the size of X " ... - "in the dimension along which variance is calculated"]); + dim = find (szx > 1, 1); + if length (dim) == 0 + dim = 1; + endif + n = szx(dim); + if (omitnan) + n = sum (! isnan (x), dim); + xn = isnan (x); + x(xn) = 0; + endif + m = sum (x, dim) ./ n; + dims = ones (1, ndx); + dims(dim) = szx(dim); + if (sprs_x) + m_exp = repmat (m, dims); + else + m_exp = m .* ones (dims); + endif + if (omitnan) + x(xn) = m_exp(xn); + endif + v = sumsq (x - m_exp, dim) ./ (n - 1 + w); + if (numel (n) == 1) + divby0 = n .* ones (size (v)) == 1; + else + divby0 = n == 1; + endif + v(divby0) = 0; + endif + + elseif (nvarg == 1) + ## Two numeric input arguments, w or weights given. + if (all_flag) + x = x(:); + if (weighted) + weights = weights(:); + wx = weights .* x; + else + weights = ones (length (x), 1); + wx = x; + endif + + if (omitnan) + xn = isnan (wx); + wx = wx(! xn); + weights = weights(! xn); + x = x(! xn); + endif + n = length (wx); + m = sum (wx) ./ sum (weights); + if (weighted) + v = sum (weights .* (abs (x - m) .^ 2)) ./ sum (weights); + else + v = sum (weights .* (abs (x - m) .^ 2)) ./ (n - 1 + w); + if (n == 1) + v = 0; + endif + endif + + else + dim = find (szx > 1, 1); + if length (dim) == 0 + dim = 1; + endif + if (! weighted) + weights = ones (szx); + wx = x; + else + if (isvector (weights)) + dims = 1:ndx; + dims([1, dim]) = [dim, 1]; + weights = zeros (szx) + permute (weights(:), dims); + endif + wx = weights .* x; + endif + n = size (wx, dim); + if (omitnan) + xn = isnan (wx); + n = sum (! xn, dim); + wx(xn) = 0; + weights(xn) = 0; + endif + m = sum (wx, dim) ./ sum (weights, dim); + dims = ones (1, ndims (wx)); + dims(dim) = size (wx, dim); + if (sprs_x) + m_exp = repmat (m, dims); + else + m_exp = m .* ones (dims); + endif + if (omitnan) + x(xn) = m_exp(xn); + endif + if (weighted) + v = sum (weights .* ((x - m_exp) .^ 2), dim) ./ sum (weights, dim); + else + v = sumsq (x - m_exp, dim) ./ (n - 1 + w); + if (numel (n) == 1) + divby0 = n .* ones (size (v)) == 1; + else + divby0 = n == 1; + endif + v(divby0) = 0; + endif + endif + + elseif (nvarg == 2) + ## Three numeric input arguments, both w or weights and dim or vecdim given. + if (vecdim_scalar_vector(1)) + if (!weighted) + weights = ones (szx); + wx = x; + else + if (isvector (weights)) + dims = 1:ndx; + dims([1, vecdim]) = [vecdim, 1]; + weights = zeros (szx) + permute (weights(:), dims); + endif + wx = weights .* x; endif - if ((dim == 1 && isrow (w)) || (dim == 2 && iscolumn (w))) - w = w.'; - elseif (dim > 2) - newdims = [ones(1, dim - 1), numel(w)]; - w = reshape (w, newdims); + n = size (wx, vecdim); + if (omitnan) + n = sum (! isnan (wx), vecdim); + xn = isnan (wx); + wx(xn) = 0; + weights(xn) = 0; + endif + m = sum (wx, vecdim) ./ sum (weights, vecdim); + dims = ones (1, ndims (wx)); + dims(vecdim) = size (wx, vecdim); + if (sprs_x) + m_exp = repmat (m, dims); + else + m_exp = m .* ones (dims); + endif + if (omitnan) + x(xn) = m_exp(xn); endif - den = sum (w); + if (weighted) + v = sum (weights .* ((x - m_exp) .^ 2), vecdim) ... + ./ sum (weights, vecdim); + else + v = sumsq (x - m_exp, vecdim); + vn = isnan (v); + v = v ./ (n - 1 + w); + if (numel (n) == 1) + divby0 = n .* ones (size (v)) == 1; + else + divby0 = n == 1; + endif + v(divby0) = 0; + v(vn) = NaN; + endif + + else + ## Weights and nonscalar vecdim specified + + ## Ignore exceeding dimensions in VECDIM + remdims = 1 : ndx; # all dimensions + vecdim(find (vecdim > ndx)) = []; + ## Calculate permutation vector + remdims(vecdim) = []; # delete dimensions specified by vecdim + nremd = numel (remdims); + + ## If all dimensions are given, it is similar to all flag + if (nremd == 0) + x = x(:); + if (weighted) + weights = weights(:); + wx = weights .* x; + else + weights = ones (length (x), 1); + wx = x; + endif - ## FIXME: Use bsxfun, rather than broadcasting, until broadcasting - ## supports diagonal and sparse matrices (Bugs #41441, #35787). - mu = sum (bsxfun (@times, w , x), dim) ./ den; - v = sum (bsxfun (@times, w, ... - bsxfun (@minus, x, mu) .^ 2), dim) ./ den; - ## mu = sum (w .* x, dim) ./ den; # automatic broadcasting - ## v = sum (w .* ((x - mu) .^ 2), dim) ./ den; + if (omitnan) + xn = isnan (wx); + wx = wx(! xn); + weights = weights(! xn); + x = x(! xn); + endif + n = length (wx); + m = sum (wx) ./ sum (weights); + if (weighted) + v = sum (weights .* (abs (x - m) .^ 2)) ./ sum (weights); + else + v = sum (weights .* (abs (x - m) .^ 2)) ./ (n - 1 + w); + if (n == 1) + v = 0; + endif + endif + + else + + ## FIXME: much of the reshaping can be skipped once octave's sum can + ## take a vecdim argument. + + ## Apply weights + if (weighted) + wx = weights .* x; + else + weights = ones (szx); + wx = x; + endif + + ## Permute to bring remaining dims forward + perm = [remdims, vecdim]; + wx = permute (wx, perm); + weights = permute (weights, perm); + x = permute (x, perm); + + ## Reshape to put all vecdims in final dimension + szwx = size (wx); + sznew = [szwx(1:nremd), prod(szwx(nremd+1:end))]; + wx = reshape (wx, sznew); + weights = reshape (weights, sznew); + x = reshape (x, sznew); + + ## Calculate var on single, squashed dimension + dim = nremd + 1; + n = size (wx, dim); + if (omitnan) + xn = isnan (wx); + n = sum (! xn, dim); + wx(xn) = 0; + weights(xn) = 0; + endif + m = sum (wx, dim) ./ sum (weights, dim); + m_exp = zeros (sznew) + m; + if (omitnan) + x(xn) = m_exp(xn); + endif + if (weighted) + v = sum (weights .* ((x - m_exp) .^ 2), dim) ./ sum (weights, dim); + else + v = sumsq (x - m_exp, dim) ./ (n - 1 + w); + if (numel (n) == 1) + divby0 = n .* ones (size (v)) == 1; + else + divby0 = n == 1; + endif + v(divby0) = 0; + endif + + ## Inverse permute back to correct dimensions + v = ipermute (v, perm); + if (nargout > 1) + m = ipermute (m, perm); + endif + endif endif endif + ## Preserve class type + if (nargout < 2) + if strcmp (outtype, "single") + v = single (v); + else + v = double (v); + endif + else + if strcmp (outtype, "single") + v = single (v); + m = single (m); + else + v = double (v); + m = double (m); + endif + endif endfunction @@ -247,6 +575,7 @@ %!assert (var (5, 99), 0) %!assert (var (5, 99, 1), 0) %!assert (var (5, 99, 2), 0) +%!assert (var ([5 3], [99 99], 2), 1) %!assert (var ([1:7], [1:7]), 3) %!assert (var ([eye(3)], [1:3]), [5/36, 2/9, 1/4], eps) %!assert (var (ones (2,2,2), [1:2], 3), [(zeros (2,2))]) @@ -254,16 +583,127 @@ %!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 single input and optional arguments "all", DIM, "omitnan") +%!test +%! x = [-10:10]; +%! y = [x;x+5;x-5]; +%! assert (var (x), 38.5); +%! assert (var (y, [], 2), [38.5; 38.5; 38.5]); +%! assert (var (y, 0, 2), [38.5; 38.5; 38.5]); +%! assert (var (y, 1, 2), ones (3,1) * 36.66666666666666, 1e-14); +%! assert (var (y, "all"), 54.19354838709678, 1e-14); +%! y(2,4) = NaN; +%! assert (var (y, "all"), NaN); +%! assert (var (y, "all", "includenan"), NaN); +%! assert (var (y, "all", "omitnan"), 55.01533580116342, 1e-14); +%! assert (var (y, 0, 2, "includenan"), [38.5; NaN; 38.5]); +%! assert (var (y, [], 2), [38.5; NaN; 38.5]); +%! assert (var (y, [], 2, "omitnan"), [38.5; 37.81842105263158; 38.5], 1e-14); + +## Tests for different weight and omitnan code paths +%!assert (var ([1 NaN 3], [1 2 3], "omitnan"), 0.75, eps) +%!assert (var ([1 2 3], [1 NaN 3], "omitnan"), 0.75, eps) +%!assert (var (magic(3), [1 NaN 3], "omitnan"), [3 12 3], eps) +%!assert (var ([1 NaN 3], [1 2 3], "omitnan", "all"), 0.75, eps) +%!assert (var ([1 NaN 3], [1 2 3], "all", "omitnan"), 0.75, eps) +%!assert (var ([1 2 3], [1 NaN 3], "omitnan", "all"), 0.75, eps) +%!assert (var ([1 NaN 3], [1 2 3], 2, "omitnan"), 0.75, eps) +%!assert (var ([1 2 3], [1 NaN 3], 2, "omitnan"), 0.75, eps) +%!assert (var (magic(3), [1 NaN 3], 1, "omitnan"), [3 12 3], eps) +%!assert (var (magic(3), [1 NaN 3], 2, "omitnan"), [0.75;3;0.75], eps) +%!assert (var ([4 4; 4 6; 6 6], [1 3], 2, 'omitnan'), [0;0.75;0], eps) +%!assert (var ([4 NaN; 4 6; 6 6], [1 2 3], 1, 'omitnan'), [1 0]) +%!assert (var ([4 NaN; 4 6; 6 6], [1 3], 2, 'omitnan'), [0;0.75;0], eps) +%!assert (var (3*reshape(1:18, [3 3 2]), [1 2 3], 1, 'omitnan'), ones(1,3,2)*5) +%!assert (var (reshape(1:18, [3 3 2]), [1 2 3], 2, 'omitnan'), 5*ones(3,1,2)) +%!assert (var (3*reshape(1:18, [3 3 2]), ones (3,3,2), [1 2], 'omitnan'), ... +%! 60 * ones(1,1,2)) +%!assert (var (3*reshape(1:18, [3 3 2]), ones (3,3,2), [1 4], 'omitnan'), ... +%! 6 * ones(1,3,2)) +%!assert (var (6*reshape(1:18, [3 3 2]), ones (3,3,2), [1:3], 'omitnan'), 969) +%!test +%! x = reshape(1:18, [3 3 2]); +%! x([2, 14]) = NaN; +%! w = ones (3,3,2); +%! assert (var (16*x, w, [1:3], 'omitnan'), 6519); +%!test +%! x = reshape(1:18, [3 3 2]); +%! w = ones (3,3,2); +%! w([2, 14]) = NaN; +%! assert (var (16*x, w, [1:3], 'omitnan'), 6519); + +## Test input case insensitivity +%!assert (var ([1 2 3], "aLl"), 1); +%!assert (var ([1 2 3], "OmitNan"), 1); +%!assert (var ([1 2 3], "IncludeNan"), 1); + +## Test dimension indexing with vecdim in n-dimensional arrays +%!test +%! x = repmat ([1:20;6:25], [5, 2, 6, 3]); +%! assert (size (var (x, 0, [3 2])), [10, 1, 1, 3]); +%! assert (size (var (x, 1, [1 2])), [1, 1, 6, 3]); +%! assert (size (var (x, [], [1 2 4])), [1, 1, 6]); +%! assert (size (var (x, 0, [1 4 3])), [1, 40]); +%! assert (size (var (x, [], [1 2 3 4])), [1, 1]); + +## Test matrix with vecdim, weighted, matrix weights, omitnan +%!assert (var (3*magic(3)), [63 144 63]) +%!assert (var (3*magic(3), 'omitnan'), [63 144 63]) +%!assert (var (3*magic(3), 1), [42 96 42]) +%!assert (var (3*magic(3), 1, 'omitnan'), [42 96 42]) +%!assert (var (3*magic(3), ones(1,3), 1), [42 96 42]) +%!assert (var (3*magic(3), ones(1,3), 1, 'omitnan'), [42 96 42]) +%!assert (var (2*magic(3), [1 1 NaN], 1, 'omitnan'), [25 16 1]) +%!assert (var (3*magic(3), ones(3,3)), [42 96 42]) +%!assert (var (3*magic(3), ones(3,3), 'omitnan'), [42 96 42]) +%!assert (var (3*magic(3), [1 1 1; 1 1 1; 1 NaN 1], 'omitnan'), [42 36 42]) +%!assert (var (3*magic(3), ones(3,3), 1), [42 96 42]) +%!assert (var (3*magic(3), ones(3,3), 1, 'omitnan'), [42 96 42]) +%!assert (var (3*magic(3), [1 1 1; 1 1 1; 1 NaN 1], 1, 'omitnan'), [42 36 42]) +%!assert (var (3*magic(3), ones(3,3), [1 4]), [42 96 42]) +%!assert (var (3*magic(3), ones(3,3), [1 4], 'omitnan'), [42 96 42]) +%!assert (var (3*magic(3), [1 1 1; 1 1 1; 1 NaN 1],[1 4],'omitnan'), [42 36 42]) + +## Test results with vecdim in n-dimensional arrays and "omitnan" +%!test +%! x = repmat ([1:20;6:25], [5, 2, 6, 3]); +%! v = repmat (33.38912133891213, [10, 1, 1, 3]); +%! assert (var (x, 0, [3, 2]), v, 1e-14); +%! v = repmat (33.250, [10, 1, 1, 3]); +%! assert (var (x, 1, [3, 2]), v, 1e-14); +%! x(2,5,6,3) = NaN; +%! v(2,1,1,3) = NaN; +%! assert (var (x, 1, [3, 2]), v, 4e-14); +%! v = repmat (33.38912133891213, [10 1 1 3]); +%! v(2,1,1,3) = NaN; +%! assert (var (x, [], [3, 2]), v, 4e-14); +%! v(2,1,1,3) = 33.40177912169048; +%! assert (var (x, [], [3, 2], "omitnan"), v, 4e-14); + +## Testing weights vector & arrays +%!assert (var (ones (2,2,2), [1:2], 3), [(zeros (2, 2))]); +%!assert (var (magic (3), [1:9], "all"), 6.666666666666667, 1e-14); + +## Test exceeding dimensions +%!assert (var (ones (2,2), [], 3), zeros (2,2)); +%!assert (var (ones (2,2,2), [], 99), zeros (2,2,2)); +%!assert (var (magic (3), [], 3), zeros (3,3)); +%!assert (var (magic (3), [], 1), [7, 16, 7]); +%!assert (var (magic (3), [], [1 3]), [7, 16, 7]); +%!assert (var (magic (3), [], [1 99]), [7, 16, 7]); + ## Test empty inputs %!assert (var ([]), NaN) +%!assert (class (var (single ([]))), "single") %!assert (var ([],[],1), NaN(1,0)) %!assert (var ([],[],2), NaN(0,1)) %!assert (var ([],[],3), []) -%!assert (var (ones (0,1)), NaN) +%!assert (class (var (single ([]), [], 1)), "single") %!assert (var (ones (1,0)), NaN) %!assert (var (ones (1,0), [], 1), NaN(1,0)) %!assert (var (ones (1,0), [], 2), NaN) %!assert (var (ones (1,0), [], 3), NaN(1,0)) +%!assert (class (var (ones (1, 0, "single"), [], 1)), "single") %!assert (var (ones (0,1)), NaN) %!assert (var (ones (0,1), [], 1), NaN) %!assert (var (ones (0,1), [], 2), NaN(0,1)) @@ -273,8 +713,11 @@ %!assert (var (ones (1,3,0,2), [], 2), NaN(1,1,0,2)) %!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 +%! [~, m] = var ([]); +%! assert (m, NaN); -## Test second output +## Test optional mean output %!test <*62395> %! [~, m] = var (13); %! assert (m, 13); @@ -289,9 +732,9 @@ %! [~, m] = var ([1, 2, 3; 3 2 1], [], 3); %! assert (m, [1 2 3; 3 2 1]); -## 2nd output, weighted inputs, vector dims +## Test mean output, weighted inputs, vector dims %!test <*62395> -%! [~, m] = var(5,99); +%! [~, m] = var (5,99); %! assert (m, 5); %! [~, m] = var ([1:7], [1:7]); %! assert (m, 5); @@ -303,19 +746,38 @@ %! assert (m, 2.5, eps); %! [~, m] = var (reshape ([1:8], 2, 2, 2), 0, [1 3]); %! assert (m, [3.5, 5.5], eps); +%!test +%! [v, m] = var (4 * eye (2), [1, 3]); +%! assert (v, [3, 3]); +%! assert (m, [1, 3]); -## 2nd output, empty inputs +## Test mean output, empty inputs, omitnan %!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)); +#%! [~, 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 mean output, nD array +%!test <*62395> +%! x = repmat ([1:20;6:25], [5, 2, 6, 3]); +%! [~, m] = var (x, 0, [3 2]); +%! assert (m, mean (x, [3 2])); +%! [~, m] = var (x, 0, [1 2]); +%! assert (m, mean (x, [1 2])); +%! [~, m] = var (x, 0, [1 3 4]); +%! assert (m, mean (x, [1 3 4])); +%!test +%! x = repmat ([1:20;6:25], [5, 2, 6, 3]); +%! x(2,5,6,3) = NaN; +%! [~, m] = var (x, 0, [3 2], "omitnan"); +%! assert (m, mean (x, [3 2], "omitnan")); ## Test Inf and NaN inputs %!test <*63203> @@ -440,27 +902,64 @@ %! [v, m] = var (sparse (4 * eye (2)), [1, 3]); %! assert (full (v), [3, 3]); %! assert (full (m), [1, 3]); - -%!test <63291> +%!test<*63291> %! [v, m] = var (sparse (eye (2))); %! assert (issparse (v)); %! assert (issparse (m)); -%!test <63291> +%!test<*63291> %! [v, m] = var (sparse (eye (2)), [1, 3]); %! assert (issparse (v)); %! assert (issparse (m)); + ## Test input validation %!error <Invalid call> var () -%!error <X must be a numeric> var (['A'; 'B']) -%!error <W must be 0> var ([1 2 3], 2) -%!error <W must be .* a vector of positive integers> var ([1 2], [-1 0]) -%!error <W must be .* a vector of positive integers> var ([1 2], eye (2)) -%!error <W must be either 0 or 1> var (ones (2, 2), [1 2], [1 2]) -%!error <W must be either 0 or 1> var ([1 2], [1 2], 'all') -%!error <the length of W must be> var ([1 2], [1 2 3]) -%!error <the length of W must be> var (1, [1 2]) -%!error <the length of W must be> var ([1 2], [1 2], 1) +%!error <Invalid call> var (1, 2, "omitnan", 3) +%!error <Invalid call> var (1, 2, 3, 4) +%!error <Invalid call> var (1, 2, 3, 4, 5) +%!error <Invalid call> var (1, "foo") +%!error <Invalid call> var (1, [], "foo") +%!error <normalization scalar must be either 0 or 1> var ([1 2 3], 2) +%!error <normalization scalar must be either 0 or 1> var ([1 2], 2, "all") +%!error <normalization scalar must be either 0 or 1> var ([1 2],0.5, "all") +%!error <weights must not contain any negative values> var (1, -1) +%!error <weights must not contain any negative values> var (1, [1 -1]) +%!error <weights must not contain any negative values> ... +%! var ([1 2 3], [1 -1 0]) +%!error <X must be a numeric vector or matrix> var ({1:5}) +%!error <X must be a numeric vector or matrix> var ("char") +%!error <X must be a numeric vector or matrix> var (['A'; 'B']) %!error <DIM must be a positive integer> var (1, [], ones (2,2)) +%!error <DIM must be a positive integer> var (1, 0, 1.5) +%!error <DIM must be a positive integer> var (1, [], 0) %!error <DIM must be a positive integer> var (1, [], 1.5) -%!error <DIM must be a positive integer> var (1, [], 0) +%!error <DIM must be a positive integer> var ([1 2 3], [], [-1 1]) +%!error <VECDIM must contain non-repeating positive integers> ... +%! var (repmat ([1:20;6:25], [5 2 6 3]), 0, [1 2 2 2]) +%!error <weight matrix or array does not match X in size> ... +%! var ([1 2], eye (2)) +%!error <weight matrix or array does not match X in size> ... +%! var ([1 2 3 4], [1 2; 3 4]) +%!error <weight matrix or array does not match X in size> ... +%! var ([1 2 3 4], [1 2; 3 4], 1) +%!error <weight matrix or array does not match X in size> ... +%! var ([1 2 3 4], [1 2; 3 4], [2 3]) +%!error <weight matrix or array does not match X in size> ... +%! var (ones (2, 2), [1 2], [1 2]) +%!error <weight matrix or array does not match X in size> ... +%! var ([1 2 3 4; 5 6 7 8], [1 2 1 2 1; 1 2 1 2 1], 1) +%!error <weight matrix or array does not match X in size> ... +%! var (repmat ([1:20;6:25], [5 2 6 3]), repmat ([1:20;6:25], [5 2 3]), [2 3]) +%!error <weight vector length does not match> var ([1 2 3; 2 3 4], [1 3 4]) +%!error <weight vector length does not match> var ([1 2], [1 2 3]) +%!error <weight vector length does not match> var (1, [1 2]) +%!error <weight vector length does not match> var ([1 2 3; 2 3 4], [1 3 4], 1) +%!error <weight vector length does not match> var ([1 2 3; 2 3 4], [1 3], 2) +%!error <weight vector length does not match> var ([1 2], [1 2], 1) +%!error <'all' flag cannot be used with DIM or VECDIM options> ... +%! var (1, [], 1, "all") +%!error <weight vector element count does not match X> ... +%! var ([1 2 3; 2 3 4], [1 3], "all") +%!error <weight matrix or array does not match X in size> ... +%! var (repmat ([1:20;6:25], [5 2 6 3]), repmat ([1:20;6:25], [5 2 3]), "all") +