Bins data on base 2 using Preston's method

    Parameters
    ----------
    data : array-like
        Data to be binned
    max_num :  float
        The maximum upper value of the data

    Returns
    -------
    tuple
        (binned_data, bin_edges)

    Notes
    -----
    Uses Preston's method of binning, which has exclusive lower boundaries and
    inclusive upper boundaries. Densities are not split between bins.

    Examples
    --------

    >>> import macroeco.compare as comp
    >>> import numpy as np

    >>> # Load some data and get Preston bins
    >>> data = np.array([1, 1, 1, 1, 4, 5, 6, 7, 12, 34, 56])
    >>> comp.preston_bin(data, np.max(data))
    (array([4, 0, 1, 3, 1, 0, 2]),
    array([  1.,   2.,   3.,   5.,   9.,  17.,  33.,  65.]))

    References
    ----------
    .. [#]
       Preston, F. (1962). The canonical distribution of commonness and rarity.
       Ecology, 43, 185-215

def preston_bin(data, max_num):
    """
    Bins data on base 2 using Preston's method

    Parameters
    ----------
    data : array-like
        Data to be binned
    max_num :  float
        The maximum upper value of the data

    Returns
    -------
    tuple
        (binned_data, bin_edges)

    Notes
    -----
    Uses Preston's method of binning, which has exclusive lower boundaries and
    inclusive upper boundaries. Densities are not split between bins.

    Examples
    --------

    >>> import macroeco.compare as comp
    >>> import numpy as np

    >>> # Load some data and get Preston bins
    >>> data = np.array([1, 1, 1, 1, 4, 5, 6, 7, 12, 34, 56])
    >>> comp.preston_bin(data, np.max(data))
    (array([4, 0, 1, 3, 1, 0, 2]),
    array([  1.,   2.,   3.,   5.,   9.,  17.,  33.,  65.]))

    References
    ----------
    .. [#]
       Preston, F. (1962). The canonical distribution of commonness and rarity.
       Ecology, 43, 185-215

    """

    log_ub = np.ceil(np.log2(max_num))

    # Make an exclusive lower bound in keeping with Preston
    if log_ub == 0:
        boundaries = np.array([0, 1])
    elif log_ub == 1:
        boundaries = np.arange(1, 4)
    else:
        boundaries = 2 ** np.arange(0, log_ub + 1)
        boundaries = np.insert(boundaries, 2, 3)
        boundaries[3:] = boundaries[3:] + 1

    hist_data = np.histogram(data, bins=boundaries)
    return hist_data