# HG changeset patch
# User Lars Winterfeld <lars@winterfeld.de>
# Date 1448040550 -3600
#      Fri Nov 20 18:29:10 2015 +0100
# Node ID 7bb3f8f0ad54ebaadb84b0696d11c98ebd7bcc89
# Parent  948f45a2931e139960a034dc750a8893d8dbdd9f
add support for the histogram and histcounts functions

more sophisticated version of histc and hist.
not yet implemented: category cell arrays as input.

diff -r 948f45a2931e -r 7bb3f8f0ad54 scripts/plot/draw/histogram.m
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/scripts/plot/draw/histogram.m	Fri Nov 20 18:29:10 2015 +0100
@@ -0,0 +1,116 @@
+## Copyright (C) 1996-2015 John W. Eaton <jwe@octave.org>
+##
+## 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 <http://www.gnu.org/licenses/>.
+
+## -*- texinfo -*-
+## @deftypefn  {Function File} {} histogram (@var{x})
+## @deftypefnx {Function File} {} histogram (@var{x}, @var{nbins})
+## @deftypefnx {Function File} {} histogram (@var{x}, @var{edges})
+## @deftypefnx {Function File} {} histogram (@dots{}, @var{name}, @var{value})
+## @deftypefnx {Function File} {} histogram (@var{hax}, @dots{})
+## @deftypefnx {Function File} {h =} histogram (@dots{})
+## Plot a histogram.
+##
+## When called with only one argument, @var{x}, plot a histogram of these
+## values. If @var{x} is not a vector, it is treated as @code{@var{x}(:)}.
+##
+## If a second scalar argument, @var{nbins}, is provided, create @var{nbins}
+## bin with uniform width.
+##
+## If a second vector argument, @var{edges}, is provided, it determines the
+## number and width of the bins. @code{@var{n}(k)} will contain the number of
+## elements in @var{x} for which
+## @code{@var{edges}(k) <= @var{x} < @var{edges}(k+1)}.
+## The last bin also includes the right edge.
+##
+## @var{name}, @var{value} pairs can be used to change the behavior of the
+## binning algorithm. Please see the documentation of @code{histcounts} for the
+## meaining of the following @var{name}s: "BinLimits", "BinMethod", "BinWidth"
+## and "Normalization". You can also provide plot options as @var{name},
+## @var{value} pairs to modify the histogram's appearance. See the examples
+## below.
+##
+## Any previous argument combinations can be preceded by an axes handle,
+## @var{hax}. In this case, plot into this axis. Otherwise, use the current
+## axes handle returned by @code{gca}.
+##
+## If an output argument, @var{h}, is requested, it currently is the used axes
+## handle.
+##
+## @example
+## @group
+## histogram (randn (1, 1000))
+## @end group
+## @end example
+##
+## @example
+## @group
+## histogram (randn (1, 1000), "BinMethod", "scott", "Normalization", "pdf")
+## hold on
+## x = linspace (-5, 5, 1000); plot (x, exp (-x.^2/2)/sqrt (2*pi), "r");
+## hold off
+## @end group
+## @end example
+##
+## @example
+## @group
+## histogram (linspace (0, 3, 1000), 23, "BinLimits", [1.0, 2.0])
+## @end group
+## @end example
+##
+## @example
+## @group
+## histogram (randn (1, 1000), "BinWidth", 0.5, "facecolor", "r", "edgecolor", "b")
+## @end group
+## @end example
+##
+## @seealso{histcounts, hist}
+## @end deftypefn
+
+## Author: Lars Winterfeld
+
+function h = histogram (varargin)
+  [hax, varargin, nargin] = __plt_get_axis_arg__ ("hist", varargin{:});
+
+  if (nargin < 1)
+    print_usage();
+  endif
+
+  if (isempty (hax))
+    hax = gca ();
+  endif
+
+  [n, edges, bin, _extraargs] = histcounts (varargin{:});
+
+  barfunc = @bar;
+  if numel(_extraargs) > 1
+    for iarg=2:numel(_extraargs)
+      if (strcmp (_extraargs{iarg-1}, "Orientation"))
+        if (strcmp (_extraargs{iarg}, "horizontal"))
+          barfunc = @barh;
+        endif
+        _extraargs([iarg-1,iarg]) = [];
+      endif
+    endfor
+  endif
+
+  hax2 = barfunc (hax, edges(1:end-1), n, "histc", _extraargs{:});
+
+  if (nargout > 0)
+    h = hax2;
+  endif
+endfunction
diff -r 948f45a2931e -r 7bb3f8f0ad54 scripts/statistics/base/histcounts.m
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/scripts/statistics/base/histcounts.m	Fri Nov 20 18:29:10 2015 +0100
@@ -0,0 +1,328 @@
+## Copyright (C) 1996-2015 John W. Eaton <jwe@octave.org>
+##
+## 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 <http://www.gnu.org/licenses/>.
+
+## -*- texinfo -*-
+## @deftypefn  {Function File} {[@var{n}, @var{edges}] =} histcounts (@var{x})
+## @deftypefnx {Function File} {[@var{n}, @var{edges}] =} histcounts (@var{x}, @var{nbins})
+## @deftypefnx {Function File} {[@var{n}, @var{edges}] =} histcounts (@var{x}, @var{edges})
+## @deftypefnx {Function File} {[@var{n}, @var{edges}, @var{bin}] =} histcounts (@dots{})
+## @deftypefnx {Function File} {[@dots{}] =} histcounts (@dots{}, @var{name}, @var{value})
+## Compute histogram counts.
+##
+## Input arguments:
+## @var{x} contains the data to be sorted into histogram bins. If @var{x} is
+## not a vector, it is treated as @code{@var{x}(:)}.
+## 
+## When called with only one argument, automatically determine the number and
+## width of bins.
+##
+## If a second scalar argument, @var{nbins}, is provided, create @var{nbins}
+## bin with uniform width.
+##
+## If a second vector argument, @var{edges}, is provided, it determines the
+## number and width of the bins. @code{@var{n}(k)} will contain the number of
+## elements in @var{x} for which
+## @code{@var{edges}(k) <= @var{x} < @var{edges}(k+1)}.
+## The last bin also includes the right edge.
+##
+## @var{name}, @var{value} pairs can be used to change the behaivor of the
+## binning algorithm.
+##
+## If a @var{name} is "BinLimits", @var{value} must be @code{[bmin, bmax]},
+## where x is truncated such that @code{@var{bmin} <= @var{x} <= @var{bmax}}.
+##
+## If a @var{name} is "BinWidth", @var{value} must be real number defining
+## the bin width for all bins.
+##
+## If a @var{name} is "BinMethod", @var{value} must be one of the following:
+## @group
+## "auto" or "fd": choose the bin width according to Freedman-Diaconis' choice
+## as @code{h = 2*iqr (@var{x}) / length (@var{x})^(1/3)}
+## @end group
+## @group
+## "integers": center bins on integers and choose bin width of 1.
+## @end group
+## @group
+## "scott": choose the bin width according to Scott's rule as
+## @code{h = 3.5*std (@var{x}) / length (@var{x})^(1/3)}
+## @end group
+## @group
+## "sqrt": choose number of bins as @code{ceil( sqrt( length (@var{x})))}
+## @end group
+## @group
+## "sturges": choose number of bins according to Sturges' rule as
+## @code{nbins = ceil (1 + log2( length (@var{x})))}
+## @end group
+##
+## If a @var{name} is "Normalization", @var{value} determines what @var{n}
+## contains and must be on of the following:
+## @group
+## "count" (default): number of elements per bin.
+## @code{sum(@var{n}) == length (@var{x})}
+## @end group
+## @group
+## "probability": (number of elements per bin) / (total number of elements).
+## @code{sum(@var{n}) == 1}
+## @end group
+## @group
+## "countdensity": (number of elements per bin) ./ (bin width).
+## @end group
+## @group
+## "pdf" (probability density function estimate):
+## (fraction of elements per bin) ./ (bin width).
+## @end group
+## @group
+## "cumcount": each @code{@var{n}(k)} contains the number of elements in the
+## @var{k}th and and all previous bins. @code{@var{n}(end) == length (@var{x})}
+## @end group
+## @group
+## "cdf" (cumulative density function estimate): each @code{@var{n}(k)} contains
+## the fraction of elements in the @var{k}th and and all previous bins.
+## @code{@var{n}(end) == 1}
+## @end group
+##
+## Output arguments:
+## @var{n} is a vector containing the number of elements in each bin.
+## @var{edges} is a vector containing the cutoffs between the bins.
+## @code{edges(1)} is the left edge of the first bin, and @code{edges(end)} is
+## the right edge of the last bin. If a third output argument, @var{bin}, is
+## requested, it is vector with the same length as @var{x} containing the index
+## of the bin that each element was assiged to.
+##
+## @example
+## @group
+## n = histcounts (randn (1, 1000))
+## @end group
+## @end example
+##
+## @example
+## @group
+## [n, edges] = histcounts (randn (1, 1000), "BinMethod", "scott", 
+## "Normalization", "pdf")
+## @end group
+## @end example
+##
+## @example
+## @group
+## [n, edges, bin] = histcounts (randn (1, 1000), [-4, 0:0.5:4])
+## @end group
+## @end example
+##
+## @example
+## @group
+## [n, edges, bin] = histcounts (linspace (0, 3, 1000), 23, "BinLimits",
+## [1.0, 2.0])
+## @end group
+## @end example
+##
+## @seealso{histogram, histc}
+## @end deftypefn
+
+## Author: Lars Winterfeld
+
+function [n, edges, bin, _extraargs] = histcounts (varargin)
+  if (nargin < 1)
+    print_usage();
+  endif
+
+  x = varargin{1};
+  x = double (x(:)); % make sure 'x' is a row vector
+  if (! (isreal (x) && isvector (x)) )
+    error ("histcounts: X must be real-valued vector");
+  endif
+  if (! all( isfinite (x)))
+    warning ("histcounts: Ignoring non-finite (Inf, -Inf, NaN) elements of X.");
+    x = x(isfinite (x));
+  endif
+
+  nbins = 0;
+  edges = [];
+  binwidths = [];
+  binwidth = 0;
+  binmethod = "auto";
+  normalization = "count";
+  _extraargs = {};
+  if (nargin > 1)
+    iarg = 2;
+    if (isscalar (varargin{iarg}) && varargin{iarg} == fix (varargin{iarg}))
+      nbins = varargin{iarg++};
+      if (nbins < 1)
+        error ("histcounts: NBINS must be at least 1");
+      endif
+    elseif (isnumeric (varargin{iarg}) && isvector (varargin{iarg})) % will be false if varargin{iarg} is empty
+      edges = varargin{iarg++};
+      if (! isreal (edges))
+        error ("histcounts: EDGES must be real-valued (not complex)");
+      endif
+      if (length (edges) <= 1)
+        error ("histcounts: EDGES must contain at least two elements");
+      endif
+      % make sure 'edges' is sorted
+      if (! issorted (edges) || edges(1) > edges(end))
+        warning ("histcounts: values were not sorted on input");
+        edges = sort (edges);
+      endif
+      edges = double (edges(:))';
+      binwidths = diff (edges);
+    endif
+
+    extraswitch = false;
+    while (iarg <= nargin)
+      switch varargin{iarg}
+        case "BinLimits"
+          binlimits = varargin{++iarg};
+          if (! isreal (binlimits) || ! isvector (binlimits))
+            error ("histcounts: BinLimits must be followed by real-valued [bmin, bmax]");
+          endif
+          x = x(binlimits(1) <= x & x <= binlimits(2));
+          if (isempty (x))
+            error ("histcounts: BinLimits left no elements");
+          endif
+        case "BinMethod"
+          binmethod = varargin{++iarg};
+          if (! ischar (binmethod))
+            error ("histcounts: BinMethod must be followed by a string");
+          endif
+        case "BinWidth"
+          binwidth = varargin{++iarg};
+          if (! isreal (binmethod))
+            error ("histcounts: BinWidth must be followed by a real number");
+          endif
+        case "Normalization"
+          normalization = varargin{++iarg};
+          if (! ischar (normalization))
+            error ("histcounts: Normalization must be followed by a string");
+          endif
+        otherwise
+          if (ischar (varargin{iarg}))
+            extraswitch = true;
+          endif
+          if (extraswitch)
+            _extraargs{numel (_extraargs)+1} = varargin{iarg};
+          else
+            error (sprintf ("argument number %i not understood", iarg));
+          endif
+      endswitch
+      ++iarg;
+    endwhile
+  endif
+
+  % if edges were not provided as argument, we will choose them to be equidistant
+  % if edges were provided as argument, we check if they are equidistant within narrow numerical tollerance
+  equidistant = isempty (edges) || max(binwidths) / min(binwidths) - 1 < 1e-12;
+
+  Nx = length (x);
+  min_val = min (x);
+  max_val = max (x);
+
+  if (binwidth <= 0 && isempty (edges))
+    switch (binmethod)
+      case {"auto", "fd"} % meaning of 'auto' not documented, we choose 'fd'
+        binwidth = 2*iqr (x) / Nx^(1/3); % Freedman-Diaconis' choice
+      case "integers"
+        edges = (floor (min_val)-0.5):(ceil (max_val)+0.5);
+      case "scott"
+        binwidth = 3.5*std (x) / Nx^(1/3); % Scott's rule
+      case "sqrt"
+        nbins = ceil( sqrt( Nx)); % sqrt(Nx) used by, e.g., MS Excel
+      case "sturges"
+        nbins = ceil (1 + log2( Nx)); % Sturges' rule
+      otherwise
+        error ("histcounts: BinMethod unknown");
+    endswitch
+  endif
+
+  % at this point, at least one of 'binwidth', 'nbins' or 'edges' is non-zero
+  if (isempty (edges))
+    if (nbins == 0)
+      nbins = round ( (max_val - min_val) / binwidth );
+      if (nbins == 0)
+        nbins = 1; % make sure we have at least one bin
+      endif
+    endif
+    edges = min_val + (max_val - min_val) * (0:nbins) / nbins;
+  elseif (nbins == 0)
+    nbins = length (edges) - 1;
+  endif
+  if isempty (binwidths) % update binwidths (since edges might have changed)
+    binwidths = diff (edges);
+  endif
+
+  if (equidistant)
+    bin = 1 + floor ((x - min_val) / binwidths(1)); % compute bin indices in O(Nx)
+  else
+    bin = lookup (edges, x); % compute bin indices in O(Nx*log(nbins)), all x < edges(1) will be zero
+  endif
+  bin(x > edges(end)) = 0;                % set too large entries to zero
+  idx = bin(bin != 0);                    % remove all out-of-scope entries
+  n = accumarray (idx, 1, [nbins+1, 1])'; % cout number of identical indices to get the histogram, in O(Nx)
+  n = [n(1:end-2), n(end-1)+n(end)];      % last bin should also contain right edge
+  % xx = edges(1:end-1) + binwidths/2;    % place bin centers half-way between the 'edges'
+
+  switch (normalization)
+    case "count"
+      ; % already done
+    case "probability"
+      n /= sum (n);
+    case "countdensity"
+      n ./= binwidths;
+    case "pdf"
+      n /= sum (n);
+      n ./= binwidths;
+    case "cumcount"
+      n = cumsum (n);
+    case "cdf"
+      n = cumsum (n);
+      n /= n(end);
+    otherwise
+      warning ("histcounts: Normalization method not understood. Falling back to 'count'.");
+  endswitch
+endfunction
+
+%!test
+%! ein = [0, 3, 6, 7, 10, 11, 12];
+%! [n, edges] = histcounts ([1, 2, 8, 9], ein);
+%! assert (n, [2 0 0 2 0 0]);
+%! assert (edges, ein);
+%!test
+%! n = histcounts ([7, 4, 3, 7, 7, 3], 3:7);
+%! assert (n, [2 1 0 3]);
+%!test
+%! [n, edges, bin] = histcounts (rand (1, 42));
+%! assert (length (bin), 42);
+%!test
+%! n = histcounts (rand (1, 42), "Normalization", "count");
+%! assert (sum(n), 42);
+%!test
+%! n = histcounts (rand (1, 42), "Normalization", "probability");
+%! assert (sum(n), 1);
+%!test
+%! n = histcounts (rand (1, 42), "Normalization", "cumcount");
+%! assert (n(end), 42);
+%!test
+%! n = histcounts (rand (1, 42), "Normalization", "cdf");
+%! assert (n(end), 1);
+
+%!error histcounts ()
+%!error histcounts ([1+i, 2, 3+3i])
+%!error histcounts (0:9, 0)
+%!error histcounts (0:9, 1.5)
+%!error histcounts (0:9, 23, 42)
+%!error histcounts (0:9, 'BinWidth', 'fail')
+%!error histcounts (6:66, 'BinLimits', [2, 3])
+%!error histcounts (0:99, 42, 'BinMethod', 'fail')
