MAINT: Cleanup for histogram bin estimator selection

Private function with superfluous checks was removed. Estimator
Estimator functions are now semi-public and selection is simplified
within histogram itself.

1 files changed, 183 insertions, 126 deletions

diff --git a/numpy/lib/function_base.py b/numpy/lib/function_base.py
index 66213c5e0..76c0e5c00 100644
--- a/numpy/lib/function_base.py
+++ b/numpy/lib/function_base.py
@@ -76,150 +76,191 @@ def iterable(y):
     return True
 
 
-def _hist_optim_numbins_estimator(a, estimator, data_range=None, data_weights=None):
+def _hist_bin_sqrt(x):
     """
-    A helper function to be called from ``histogram`` to deal with
-    estimating optimal number of bins.
+    Square root histogram bin estimator.
 
-    A description of the estimators can be found at
-    https://en.wikipedia.org/wiki/Histogram#Number_of_bins_and_width
+    Used by many programs for its simplicity.
 
     Parameters
     ----------
-    a : array_like
-        The data with which to estimate the number of bins
-    estimator: str
-        If ``estimator`` is one of ['auto', 'fd', 'scott', 'doane',
-        'rice', 'sturges', 'sqrt'], this function will choose the
-        appropriate estimation method and return the optimal number of
-        bins it calculates.
-    data_range: tuple (min, max)
-        The range that the data to be binned should be restricted to.
-    data_weights:
-        weights are not supported, so this field must be empty or None.
+    x : array_like
+        Input data that is to be histogrammed.
+
+    Returns
+    -------
+    n : An estimate of the optimal bin count for the given data.
     """
-    if a.size == 0:
-        return 1
+    return int(np.ceil(np.sqrt(x.size)))
 
-    if data_weights is not None:
-        raise TypeError("Automated estimation of the number of "
-                        "bins is not supported for weighted data")
 
-    if data_range is not None:
-        mn, mx = data_range
-        keep = (a >= mn)
-        keep &= (a <= mx)
-        if not np.logical_and.reduce(keep):
-            a = a[keep]
+def _hist_bin_sturges(x):
+    """
+    Sturges histogram bin estimator.
 
-    def sqrt(x):
-        """
-        Square Root Estimator
+    A very simplistic estimator based on the assumption of normality of
+    the data. This estimator has poor performance for non-normal data,
+    which becomes especially obvious for large data sets. The estimate
+    depends only on size of the data.
 
-        Used by many programs for its simplicity.
-        """
-        return np.ceil(np.sqrt(x.size))
+    Parameters
+    ----------
+    x : array_like
+        Input data that is to be histogrammed.
 
-    def sturges(x):
-        """
-        Sturges Estimator
+    Returns
+    -------
+    n : An estimate of the optimal bin count for the given data.
+    """
+    return int(np.ceil(np.log2(x.size))) + 1
 
-        A very simplistic estimator based on the assumption of normality
-        of the data. Poor performance for non-normal data, especially
-        obvious for large ``x``. Depends only on size of the data.
-        """
-        return np.ceil(np.log2(x.size)) + 1
 
-    def rice(x):
-        """
-        Rice Estimator
+def _hist_bin_rice(x):
+    """
+    Rice histogram bin estimator.
 
-        Another simple estimator, with no normality assumption. It has
-        better performance for large data, but tends to overestimate
-        number of bins. The number of bins is proportional to the cube
-        root of data size (asymptotically optimal). Depends only on size
-        of the data.
-        """
-        return np.ceil(2 * x.size ** (1.0 / 3))
+    Another simple estimator with no normality assumption. It has better
+    performance for large data than Sturges, but tends to overestimate
+    the number of bins. The number of bins is proportional to the cube
+    root of data size (asymptotically optimal). The estimate depends
+    only on size of the data.
 
-    def scott(x):
-        """
-        Scott Estimator
+    Parameters
+    ----------
+    x : array_like
+        Input data that is to be histogrammed.
 
-        The binwidth is proportional to the standard deviation of the
-        data and inversely proportional to the cube root of data size
-        (asymptotically optimal).
-        """
-        h = (24 * np.pi**0.5 / x.size)**(1.0 / 3) * np.std(x)
-        if h > 0:
-            return np.ceil(x.ptp() / h)
-        return 1
+    Returns
+    -------
+    n : An estimate of the optimal bin count for the given data.
+    """
+    return int(np.ceil(2 * x.size ** (1.0 / 3)))
 
-    def doane(x):
-        """
-        Doane's Estimator
 
-        Improved version of Sturges' formula which works better for
-        non-normal data. See
-        http://stats.stackexchange.com/questions/55134/doanes-formula-for-histogram-binning
-        """
-        if x.size > 2:
-            sg1 = np.sqrt(6.0 * (x.size - 2) / ((x.size + 1.0) * (x.size + 3)))
-            sigma = np.std(x)
-            if sigma > 0:
-                # These three operations add up to
-                # g1 = np.mean(((x - np.mean(x)) / sigma)**3)
-                # but use only one temp array instead of three
-                temp = x - np.mean(x)
-                np.true_divide(temp, sigma, temp)
-                np.power(temp, 3, temp)
-                g1 = np.mean(temp)
-                return np.ceil(1.0 + np.log2(x.size) +
-                                     np.log2(1.0 + np.absolute(g1) / sg1))
-        return 1
-
-    def fd(x):
-        """
-        Freedman Diaconis Estimator
+def _hist_bin_scott(x):
+    """
+    Scott histogram bin estimator.
 
-        The interquartile range (IQR) is used for binwidth, making this
-        variation of the Scott rule more robust, as the IQR is less
-        affected by outliers than the standard deviation. However, the
-        IQR depends on fewer points than the standard deviation, so it
-        is less accurate, especially for long tailed distributions.
+    The binwidth is proportional to the standard deviation of the data
+    and inversely proportional to the cube root of data size
+    (asymptotically optimal).
 
-        If the IQR is 0, we return 1 for the number of bins. Binwidth is
-        inversely proportional to the cube root of data size
-        (asymptotically optimal).
-        """
-        iqr = np.subtract(*np.percentile(x, [75, 25]))
+    Parameters
+    ----------
+    x : array_like
+        Input data that is to be histogrammed.
 
-        if iqr > 0:
-            h = (2 * iqr * x.size ** (-1.0 / 3))
-            return np.ceil(x.ptp() / h)
+    Returns
+    -------
+    n : An estimate of the optimal bin count for the given data.
+    """
+    h = (24 * np.pi**0.5 / x.size)**(1.0 / 3) * np.std(x)
+    if h > 0:
+        return int(np.ceil(x.ptp() / h))
+    return 1
 
-        # If iqr is 0, default number of bins is 1
-        return 1
 
-    def auto(x):
-        """
-        The FD estimator is usually the most robust method, but it tends
-        to be too small for small ``x``. The Sturges estimator is quite
-        good for small (<1000) datasets and is the default in R. This
-        method gives good off-the-shelf behaviour.
-        """
-        return max(fd(x), sturges(x))
+def _hist_bin_doane(x):
+    """
+    Doane's histogram bin estimator.
 
-    optimal_numbins_methods = {'sqrt': sqrt, 'sturges': sturges,
-                               'rice': rice, 'scott': scott, 'doane': doane,
-                               'fd': fd, 'auto': auto}
-    try:
-        estimator_func = optimal_numbins_methods[estimator.lower()]
-    except KeyError:
-        raise ValueError("{0} not a valid method for `bins`".format(estimator))
-    else:
-        # these methods return floats, np.histogram requires an int
-        return int(estimator_func(a))
+    Improved version of Sturges' formula which works better for
+    non-normal data. See
+    http://stats.stackexchange.com/questions/55134/doanes-formula-for-histogram-binning
+
+    Parameters
+    ----------
+    x : array_like
+        Input data that is to be histogrammed.
+
+    Returns
+    -------
+    n : An estimate of the optimal bin count for the given data.
+    """
+    if x.size > 2:
+        sg1 = np.sqrt(6.0 * (x.size - 2) / ((x.size + 1.0) * (x.size + 3)))
+        sigma = np.std(x)
+        if sigma > 0:
+            # These three operations add up to
+            # g1 = np.mean(((x - np.mean(x)) / sigma)**3)
+            # but use only one temp array instead of three
+            temp = x - np.mean(x)
+            np.true_divide(temp, sigma, temp)
+            np.power(temp, 3, temp)
+            g1 = np.mean(temp)
+            return int(np.ceil(1.0 + np.log2(x.size) +
+                                     np.log2(1.0 + np.absolute(g1) / sg1)))
+    return 1
+
+
+def _hist_bin_fd(x):
+    """
+    The Freedman-Diaconis histogram bin estimator.
+
+    The Freedman-Diaconis rule uses interquartile range (IQR)
+    binwidth. It is considered a variation of the Scott rule with more
+    robustness as the IQR is less affected by outliers than the standard
+    deviation. However, the IQR depends on fewer points than the
+    standard deviation, so it is less accurate, especially for long
+    tailed distributions.
+
+    If the IQR is 0, this function returns 1 for the number of bins.
+    Binwidth is inversely proportional to the cube root of data size
+    (asymptotically optimal).
+
+    Parameters
+    ----------
+    x : array_like
+        Input data that is to be histogrammed.
+
+    Returns
+    -------
+    n : An estimate of the optimal bin count for the given data.
+    """
+    iqr = np.subtract(*np.percentile(x, [75, 25]))
+
+    if iqr > 0:
+        h = (2 * iqr * x.size ** (-1.0 / 3))
+        return int(np.ceil(x.ptp() / h))
+
+    # If iqr is 0, default number of bins is 1
+    return 1
+
+
+def _hist_bin_auto(x):
+    """
+    Histogram bin estimator that uses the maximum of the
+    Freedman-Diaconis and Sturges estimators.
+
+    The FD estimator is usually the most robust method, but it tends to
+    be too small for small `x`. The Sturges estimator is quite good for
+    small (<1000) datasets and is the default in the R language. This
+    method gives good off the shelf behaviour.
+
+    Parameters
+    ----------
+    x : array_like
+        Input data that is to be histogrammed.
+
+    Returns
+    -------
+    n : An estimate of the optimal bin count for the given data.
+
+    See Also
+    --------
+    _hist_bin_fd, _hist_bin_sturges
+    """
+    return max(_hist_bin_fd(x), _hist_bin_sturges(x))
+
+
+# Private dict initialized at module load time
+_hist_bin_selectors = {'auto': _hist_bin_auto,
+                       'doane': _hist_bin_doane,
+                       'fd': _hist_bin_fd,
+                       'rice': _hist_bin_rice,
+                       'scott': _hist_bin_scott,
+                       'sqrt': _hist_bin_sqrt,
+                       'sturges': _hist_bin_sturges}
 
 
 def histogram(a, bins=10, range=None, normed=False, weights=None,
@@ -241,9 +282,9 @@ def histogram(a, bins=10, range=None, normed=False, weights=None,
 
         If `bins` is a string from the list below, `histogram` will use
         the method chosen to calculate the optimal number of bins (see
-        Notes for more detail on the estimators). For visualisation, we
-        suggest using the 'auto' option. Weighted data is not supported
-        for automated bin size selection.
+        `Notes` for more detail on the estimators). For visualisation,
+        using the 'auto' option is suggested. Weighted data is not
+        supported for automated bin size selection.
 
         'auto'
             Maximum of the 'sturges' and 'fd' estimators. Provides good
@@ -342,7 +383,7 @@ def histogram(a, bins=10, range=None, normed=False, weights=None,
         value will usually be chosen, while larger datasets will usually
         default to FD.  Avoids the overly conservative behaviour of FD
         and Sturges for small and large datasets respectively.
-        Switchover point usually happens when ``x.size`` is around 1000.
+        Switchover point is usually :math:`a.size \approx 1000`.
 
     'FD' (Freedman Diaconis Estimator)
         .. math:: h = 2 \frac{IQR}{n^{1/3}}
@@ -444,11 +485,27 @@ def histogram(a, bins=10, range=None, normed=False, weights=None,
             raise ValueError(
                 'range parameter must be finite.')
 
-
     if isinstance(bins, basestring):
-        bins = _hist_optim_numbins_estimator(a, bins, range, weights)
         # if `bins` is a string for an automatic method,
         # this will replace it with the number of bins calculated
+        if bins not in _hist_bin_selectors:
+            raise ValueError("{0} not a valid estimator for `bins`".format(bins))
+        if weights is not None:
+            raise TypeError("Automated estimation of the number of "
+                            "bins is not supported for weighted data")
+        # Make a reference to `a`
+        b = a
+        # Update the reference if the range needs truncation
+        if range is not None:
+            mn, mx = range
+            keep = (a >= mn)
+            keep &= (a <= mx)
+            if not np.logical_and.reduce(keep):
+                b = a[keep]
+        if b.size == 0:
+            bins = 1
+        else:
+            bins = _hist_bin_selectors[bins](b)
 
     # Histogram is an integer or a float array depending on the weights.
     if weights is None: