octave_packages/m/statistics/base/histc.m

   1 ## Copyright (C) 2009-2012 Søren Hauberg
   2 ## Copyright (C) 2009 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} {@var{n} =} histc (@var{x}, @var{edges})
  22 ## @deftypefnx {Function File} {@var{n} =} histc (@var{x}, @var{edges}, @var{dim})
  23 ## @deftypefnx {Function File} {[@var{n}, @var{idx}] =} histc (@dots{})
  24 ## Produce histogram counts.
  25 ##
  26 ## When @var{x} is a vector, the function counts the number of elements of
  27 ## @var{x} that fall in the histogram bins defined by @var{edges}.  This must be
  28 ## a vector of monotonically increasing values that define the edges of the
  29 ## histogram bins.  @code{@var{n}(k)} contains the number of elements in
  30 ## @var{x} for which @code{@var{edges}(k) <= @var{x} < @var{edges}(k+1)}.
  31 ## The final element of @var{n} contains the number of elements of @var{x}
  32 ## exactly equal to the last element of @var{edges}.
  33 ##
  34 ## When @var{x} is an @math{N}-dimensional array, the computation is
  35 ## carried out along dimension @var{dim}.  If not specified @var{dim} defaults
  36 ## to the first non-singleton dimension.
  37 ##
  38 ## When a second output argument is requested an index matrix is also returned.
  39 ## The @var{idx} matrix has the same size as @var{x}.  Each element of @var{idx}
  40 ## contains the index of the histogram bin in which the corresponding element
  41 ## of @var{x} was counted.
  42 ## @seealso{hist}
  43 ## @end deftypefn
  44
  45 function [n, idx] = histc (x, edges, dim)
  46
  47   if (nargin < 2 || nargin > 3)
  48     print_usage ();
  49   endif
  50
  51   if (!isreal (x))
  52     error ("histc: X argument must be real-valued, not complex");
  53   endif
  54
  55   num_edges = numel (edges);
  56   if (num_edges == 0)
  57     error ("histc: EDGES must not be empty");
  58   endif
  59
  60   if (!isreal (edges))
  61     error ("histc: EDGES must be real-valued, not complex");
  62   else
  63     ## Make sure 'edges' is sorted
  64     edges = edges(:);
  65     if (!issorted (edges) || edges(1) > edges(end))
  66       warning ("histc: edge values not sorted on input");
  67       edges = sort (edges);
  68     endif
  69   endif
  70
  71   nd = ndims (x);
  72   sz = size (x);
  73   if (nargin < 3)
  74     ## Find the first non-singleton dimension.
  75     (dim = find (sz > 1, 1)) || (dim = 1);
  76   else
  77     if (!(isscalar (dim) && dim == fix (dim))
  78         || !(1 <= dim && dim <= nd))
  79       error ("histc: DIM must be an integer and a valid dimension");
  80     endif
  81   endif
  82
  83   nsz = sz;
  84   nsz(dim) = num_edges;
  85
  86   ## the splitting point is 3 bins
  87
  88   if (num_edges <= 3)
  89
  90     ## This is the O(M*N) algorithm.
  91
  92     ## Allocate the histogram
  93     n = zeros (nsz);
  94
  95     ## Allocate 'idx'
  96     if (nargout > 1)
  97       idx = zeros (sz);
  98     endif
  99
 100     ## Prepare indices
 101     idx1 = cell (1, dim-1);
 102     for k = 1:length (idx1)
 103       idx1{k} = 1:sz(k);
 104     endfor
 105     idx2 = cell (length (sz) - dim);
 106     for k = 1:length (idx2)
 107       idx2{k} = 1:sz(k+dim);
 108     endfor
 109
 110     ## Compute the histograms
 111     for k = 1:num_edges-1
 112       b = (edges (k) <= x & x < edges (k+1));
 113       n(idx1{:}, k, idx2{:}) = sum (b, dim);
 114       if (nargout > 1)
 115         idx(b) = k;
 116       endif
 117     endfor
 118     b = (x == edges (end));
 119     n(idx1{:}, num_edges, idx2{:}) = sum (b, dim);
 120     if (nargout > 1)
 121       idx(b) = num_edges;
 122     endif
 123
 124   else
 125
 126     ## This is the O(M*log(N) + N) algorithm.
 127
 128     ## Look-up indices.
 129     idx = lookup (edges, x);
 130     ## Zero invalid ones (including NaNs). x < edges(1) are already zero.
 131     idx(! (x <= edges(end))) = 0;
 132
 133     iidx = idx;
 134
 135     ## In case of matrix input, we adjust the indices.
 136     if (! isvector (x))
 137       nl = prod (sz(1:dim-1));
 138       nn = sz(dim);
 139       nu = prod (sz(dim+1:end));
 140       if (nl != 1)
 141         iidx = (iidx-1) * nl;
 142         iidx += reshape (kron (ones (1, nn*nu), 1:nl), sz);
 143       endif
 144       if (nu != 1)
 145         ne =length (edges);
 146         iidx += reshape (kron (nl*ne*(0:nu-1), ones (1, nl*nn)), sz);
 147       endif
 148     endif
 149
 150     ## Select valid elements.
 151     iidx = iidx(idx != 0);
 152
 153     ## Call accumarray to sum the indexed elements.
 154     n = accumarray (iidx(:), 1, nsz);
 155
 156   endif
 157
 158 endfunction
 159
 160
 161 %!test
 162 %! x = linspace (0, 10, 1001);
 163 %! n = histc (x, 0:10);
 164 %! assert (n, [repmat(100, 1, 10), 1]);
 165
 166 %!test
 167 %! x = repmat (linspace (0, 10, 1001), [2, 1, 3]);
 168 %! n = histc (x, 0:10, 2);
 169 %! assert (n, repmat ([repmat(100, 1, 10), 1], [2, 1, 3]));
 170
 171 %!error histc ();
 172 %!error histc (1);
 173 %!error histc (1, 2, 3, 4);
 174 %!error histc ([1:10 1+i], 2);
 175 %!error histc (1:10, []);
 176 %!error histc (1, 1, 3);