######################################################################## ## ## Copyright (C) 1996-2023 The Octave Project Developers ## ## See the file COPYRIGHT.md in the top-level directory of this ## distribution or . ## ## This file is part of Octave. ## ## Octave is free software: you can redistribute it and/or modify it ## under the terms of the GNU General Public License as published by ## the Free Software Foundation, either version 3 of the License, or ## (at your option) any later version. ## ## Octave is distributed in the hope that it will be useful, but ## WITHOUT ANY WARRANTY; without even the implied warranty of ## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the ## GNU General Public License for more details. ## ## You should have received a copy of the GNU General Public License ## along with Octave; see the file COPYING. If not, see ## . ## ######################################################################## ## -*- texinfo -*- ## @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}. ## ## 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 ## ## @example ## mean (@var{x}) = SUM_i @var{x}(i) / N ## @end example ## ## @noindent ## where @math{N} is the number of elements in @var{x}. ## ## @end ifnottex ## ## 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}. ## ## 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. ## ## 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. ## @var{outtype} can take the following values: ## ## @table @asis ## @item @qcode{'default'} : Output is of type double, unless the input is ## single in which case the output is of type single. ## ## @item @qcode{'double'} : Output is of type double. ## ## @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 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, movmean} ## @end deftypefn function m = mean (x, varargin) if (nargin < 1 || nargin > 4) print_usage (); endif ## 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 ## 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 "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 endif 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) || ischar (x))) error ("mean: X must be either a numeric, boolean, or character array."); endif ## Process special cases for in/out size if (nvarg == 0) ## Single numeric input argument, no dimensions given. if (all_flag) x = x(:); if (omitnan) 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 else ## Find the first non-singleton dimension. (dim = find (szx != 1, 1)) || (dim = 1); n = szx(dim); if (omitnan) idx = isnan (x); n = sum (! idx, dim); x(idx) = 0; endif 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! 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 (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 if (isscalar (vecdim)) if (vecdim > ndx) m = x; else n = szx(vecdim); if (omitnan) nanx = isnan (x); n = sum (! nanx, vecdim); x(nanx) = 0; endif if (any (isa (x, {"int64", "uint64"}))) m = int64_mean (x, vecdim, n, outtype); else m = sum (x, vecdim) ./ n; endif endif 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 = 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 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); ## Calculate mean on dim1 if (omitnan) nanx = isnan (x); x(nanx) = 0; n = sum (! nanx, 1); else n = szx(1); endif 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 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; %! y = x'; %! z = [y, y+10]; %! assert (mean (x), 0); %! assert (mean (y), 0); %! assert (mean (z), [0, 10]); %!assert (mean (magic (3), 1), [5, 5, 5]) %!assert (mean (magic (3), 2), [5; 5; 5]) %!assert (mean (logical ([1 0 1 1])), 0.75) %!assert (mean (single ([1 0 1 1])), single (0.75)) %!assert (mean ([1 2], 3), [1 2]) #### 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)); %!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, "double"), out); %! assert (mean (in, "native"), uint8 (out_u8)); %! assert (class (mean (in, "native")), "uint8"); %!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]; %! assert (mean (x), 0); %! assert (mean (y, 2), [0, 5, -5]'); %! assert (mean (y, "all"), 0); %! y(2,4) = NaN; %! 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); %! m = [20 25.35 15]; %! assert (mean (z', "omitnan"), m); %! assert (mean (z, 2, "omitnan"), m'); %! assert (mean (z, 2, "native", "omitnan"), m'); %! assert (mean (z, 2, "omitnan", "native"), m'); # Test boolean input %!test %! assert (mean (true, "all"), 1); %! assert (mean (false), 0); %! assert (mean ([true false true]), 2/3, 4e-14); %! assert (mean ([true false true], 1), [1 0 1]); %! assert (mean ([true false NaN], 1), [1 0 NaN]); %! assert (mean ([true false NaN], 2), NaN); %! assert (mean ([true false NaN], 2, "omitnan"), 0.5); %! assert (mean ([true false NaN], 2, "omitnan", "native"), 0.5); ## 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]); %! assert (size (mean (x, [1 2])), [1 1 6 3]); %! assert (size (mean (x, [1 2 4])), [1 1 6]); %! assert (size (mean (x, [1 4 3])), [1 40]); %! assert (size (mean (x, [1 2 3 4])), [1 1]); ## 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,1,1,3) = NaN; %! assert (mean (x, [3 2]), m, 4e-14); %! 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 mean () %!error mean (1, 2, 3) %!error mean (1, 2, 3, 4) %!error mean (1, "all", 3) %!error mean (1, "b") %!error mean (1, 1, "foo") %!error mean ("abc", "native") %!error mean ({1:5}) %!error mean (1, ones (2,2)) %!error mean (1, 1.5) %!error mean (1, 0) %!error mean (1, []) %!error mean (repmat ([1:20;6:25], [5 2]), -1) %!error mean (repmat ([1:5;5:9], [5 2]), [1 -1]) %!error mean (1, ones(1,0)) %!error mean (1, [2 2])