######################################################################## ## ## Copyright (C) 1995-2021 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 (@var{x}) ## @deftypefnx {} {} var (@var{x}, @var{w}) ## @deftypefnx {} {} var (@var{x}, @var{w}, @var{dim}) ## @deftypefnx {} {} var (@var{x}, @var{w}, @qcode{"ALL"}) ## 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} ## $$ ## 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 ## @end group ## @end example ## ## @noindent ## where @math{N} is the length of the @var{x} vector. ## ## @end ifnottex ## If @var{x} is a matrix, compute the variance for each column and return ## them in a row vector. ## ## The argument @var{w} determines the weighting scheme to use. ## Valid values are ## ## @table @asis ## @item 0: ## normalize with @math{N-1}, provides the square root of the best unbiased ## estimator of the variance [default] ## ## @item 1: ## normalize with @math{N}, this provides the square root of the second ## moment around the mean ## ## @item a vector: ## compute the weighted variance with nonnegative scalar weights, ## the length of @var{w} must be equal to the size of @var{x} along dimension ## @var{dim} ## @end table ## ## If @math{N} is equal to 1 the value of @var{opt} is ignored and ## normalization by @math{N} is used. ## ## The optional variable @var{dim} can be used to force @code{var} to operate ## over the specified dimension. @var{dim} can either be a scalar dimension or ## a vector of non-repeating dimensions over which to operate. Dimensions are ## positive integers. When @var{dim} is a vector, @var{w} must be either 0 or ## 1, and the variance is calculated over the array slice defined by @var{dim}. ## ## Specifying dimension @qcode{"ALL"} will force @code{var} to operate ## on all elements of @var{x}, when @var{w} is either 0 or 1, and is equivalent ## to @code{var (@var{x}(:))}. ## @seealso{cov, std, skewness, kurtosis, moment} ## @end deftypefn function retval = var (x, w = 0, dim) if (nargin < 1) print_usage (); elseif (nargin < 3) dim = []; endif if (! (isnumeric (x) || islogical (x))) error ("var: X must be a numeric vector or matrix"); endif if (isempty (w)) w = 0; elseif (! isvector (w) || ! isnumeric (w) || (isvector (w) && any (w < 0)) || (isscalar (w) && (w != 0 && w != 1))) error ("var: W must be 0 or 1 or a vector of positive integers"); endif nd = ndims (x); sz = size (x); emptydimflag = false; if (isempty (dim)) emptydimflag = true; ## Compatibliity hack for empty x, ndims==2 ## Find the first non-singleton dimension. (dim = find (sz != 1, 1)) || (dim = 1); else if (! (isscalar (dim) && dim == fix (dim) && dim > 0)) if (isvector (dim) && isnumeric (dim) && all (dim > 0) && all (rem (dim, 1) == 0)) if (dim != unique (dim, "stable")) error (["var: vector DIM must contain non-repeating positive"... "integers"]); endif ## Check W if (! isscalar (w)) error ("var: W must be either 0 or 1 when DIM is a vector"); endif ## Reshape X to compute the variance over an array slice if (iscolumn (dim)) dim = transpose (dim); endif collapsed_dims = dim; dim = dim(end); ## Permute X to cluster the dimensions to collapse highest_dim = max ([nd, collapsed_dims]); perm_start = perm_end = [1:highest_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); ## Collapse the given dimensions newshape = ones (1, highest_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) && strcmp (tolower (dim), "all")) ## Check W if (! isscalar (w)) error ("var: W must be either 0 or 1 when using 'ALL' as dimension"); endif ## "ALL" equals 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'"); endif endif endif if (isempty (x)) if (emptydimflag && isequal (sz, [0 0])) retval = NaN; else output_size = sz; output_size(dim) = 1; retval = NaN(output_size); endif else n = size (x, dim); if (n == 1) if (isa (x, "single")) retval = zeros (sz, "single"); else retval = zeros (sz); endif else if (isscalar (w)) retval = sumsq (center (x, dim), dim) / (n - 1 + w); else ## Weighted variance if (length (w) != n) error (["var: the length of W must be equal to the size of X "... "along which the variance is calculated"]); else if ((dim == 1 && rows (w) == 1) || (dim == 2 && columns (w) == 1)) w = transpose (w); elseif (dim > 2) newdims = [(ones (1, (dim - 1))), (length (w))]; w = reshape (w, newdims); endif den = sum (w); mu = sum (w .* x, dim) ./ sum (w); retval = sum (w .* ((x .- mu) .^ 2), dim) / den; endif endif endif endif endfunction ##%!assert (var ([]), NaN) %!assert (var (13), 0) %!assert (var (single (13)), single (0)) %!assert (var ([1,2,3]), 1) %!assert (var ([1,2,3], 1), 2/3, eps) %!assert (var ([1,2,3], [], 1), [0,0,0]) %!assert (var ([1,2,3], [], 3), [0,0,0]) %!assert (var ([1:7], [1:7]), 3) %!assert (var ([(eye (3))], [1:3]), [.1389 .2222 .25], 1e-4) %!assert (var (ones (2,2,2), [1:2], 3), [(zeros (2,2))]) %!assert (var ([1 2; 3 4], 0, 'all'), var ([1:4])) %!assert (var (reshape ([1:8], 2, 2, 2), 0, [1 3]), [5.6667 5.6667], 1e-4) ##Test empty inputs %!assert (var ([]), NaN) %!assert (var ([],[],1), NaN(1,0)) %!assert (var ([],[],2), NaN(0,1)) %!assert (var ([],[],3), []) %!assert (var (ones (0,1)), NaN) %!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 (var (ones (0,1)), NaN) %!assert (var (ones (0,1), [], 1), NaN) %!assert (var (ones (0,1), [], 2), NaN(0,1)) %!assert (var (ones (0,1), [], 3), NaN(0,1)) %!assert (var (ones (1,3,0,2)), NaN(1,1,0,2)) %!assert (var (ones (1,3,0,2), [], 1), NaN(1,3,0,2)) %!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 input validation %!error var () %!error var (['A'; 'B']) %!error var (1, 2) %!error var ([1 2], [-1 0]) %!error var ([1 2], eye (2)) %!error var (ones (2, 2), [1 2], [1 2]) %!error var ([1 2], [1 2], 'all') %!error var ([1 2], [1 2 3]) %!error var (1, [], ones (2,2)) %!error var (1, [], 1.5) %!error var (1, [], 0)