octave_packages/m/statistics/base/median.m

   1 ## Copyright (C) 1996-2012 John W. Eaton
   2 ## Copyright (C) 2009-2010 VZLU Prague
   3 ##
   4 ## This file is part of Octave.
   5 ##
   6 ## Octave is free software; you can redistribute it and/or modify it
   7 ## under the terms of the GNU General Public License as published by
   8 ## the Free Software Foundation; either version 3 of the License, or (at
   9 ## your option) any later version.
  10 ##
  11 ## Octave is distributed in the hope that it will be useful, but
  12 ## WITHOUT ANY WARRANTY; without even the implied warranty of
  13 ## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  14 ## General Public License for more details.
  15 ##
  16 ## You should have received a copy of the GNU General Public License
  17 ## along with Octave; see the file COPYING.  If not, see
  18 ## <http://www.gnu.org/licenses/>.
  19
  20 ## -*- texinfo -*-
  21 ## @deftypefn  {Function File} {} median (@var{x})
  22 ## @deftypefnx {Function File} {} median (@var{x}, @var{dim})
  23 ## Compute the median value of the elements of the vector @var{x}.
  24 ## If the elements of @var{x} are sorted, the median is defined
  25 ## as
  26 ## @tex
  27 ## $$
  28 ## {\rm median} (x) =
  29 ##   \cases{x(\lceil N/2\rceil), & $N$ odd;\cr
  30 ##           (x(N/2)+x(N/2+1))/2, & $N$ even.}
  31 ## $$
  32 ## @end tex
  33 ## @ifnottex
  34 ##
  35 ## @example
  36 ## @group
  37 ##               x(ceil(N/2))             N odd
  38 ## median (x) =
  39 ##              (x(N/2) + x((N/2)+1))/2   N even
  40 ## @end group
  41 ## @end example
  42 ##
  43 ## @end ifnottex
  44 ## If @var{x} is a matrix, compute the median value for each
  45 ## column and return them in a row vector.  If the optional @var{dim}
  46 ## argument is given, operate along this dimension.
  47 ## @seealso{mean, mode}
  48 ## @end deftypefn
  49
  50 ## Author: jwe
  51
  52 function retval = median (x, dim)
  53
  54   if (nargin != 1 && nargin != 2)
  55     print_usage ();
  56   endif
  57
  58   if (! (isnumeric (x) || islogical (x)))
  59     error ("median: X must be a numeric vector or matrix");
  60   endif
  61
  62   if (isempty (x))
  63     error ("median: X cannot be an empty matrix");
  64   endif
  65
  66   nd = ndims (x);
  67   sz = size (x);
  68   if (nargin < 2)
  69     ## Find the first non-singleton dimension.
  70     (dim = find (sz > 1, 1)) || (dim = 1);
  71   else
  72     if (!(isscalar (dim) && dim == fix (dim))
  73         || !(1 <= dim && dim <= nd))
  74       error ("median: DIM must be an integer and a valid dimension");
  75     endif
  76   endif
  77
  78   n = sz(dim);
  79   k = floor ((n+1) / 2);
  80   if (mod (n, 2) == 1)
  81     retval = nth_element (x, k, dim);
  82   else
  83     retval = mean (nth_element (x, k:k+1, dim), dim);
  84   endif
  85   ## Inject NaNs where needed, to be consistent with Matlab.
  86   retval(any (isnan (x), dim)) = NaN;
  87
  88 endfunction
  89
  90
  91 %!test
  92 %! x = [1, 2, 3, 4, 5, 6];
  93 %! x2 = x';
  94 %! y = [1, 2, 3, 4, 5, 6, 7];
  95 %! y2 = y';
  96 %!
  97 %! assert(median (x) == median (x2) && median (x) == 3.5);
  98 %! assert(median (y) == median (y2) && median (y) == 4);
  99 %! assert(median ([x2, 2*x2]) == [3.5, 7]);
 100 %! assert(median ([y2, 3*y2]) == [4, 12]);
 101
 102 %% Test multidimensional arrays (bug #35679)
 103 %!shared a, b, x, y
 104 %! rand ("seed", 2);
 105 %! a = rand (2,3,4,5);
 106 %! b = rand (3,4,6,5);
 107 %! x = sort (a, 4);
 108 %! y = sort (b, 3);
 109 %!assert (median (a, 4), x(:, :, :, 3));
 110 %!assert (median (b, 3), (y(:, :, 3, :) + y(:, :, 4, :))/2);
 111
 112 %!assert(median (single([1,2,3])), single(2));
 113 %!assert(median ([1,2,NaN;4,5,6;NaN,8,9]), [NaN, 5, NaN]);
 114
 115 %% Test input validation
 116 %!error median ();
 117 %!error median (1, 2, 3);
 118 %!error median ({1:5});
 119 %!error median (['A'; 'B']);
 120 %!error median (1, ones(2,2));
 121 %!error median (1, 1.5);
 122 %!error median (1, 0);
 123