Returns d-1 dimensional histogram containing percentile of values along axis
        if inclusive=True, will report bin center of first bin for which percentile% of data lies in or below the bin
                    =False, ... data lies strictly below the bin
        10% percentile is calculated as: value at least 10% data is LOWER than

def percentile(self, percentile, axis, inclusive=True):
        """Returns d-1 dimensional histogram containing percentile of values along axis
        if inclusive=True, will report bin center of first bin for which percentile% of data lies in or below the bin
                    =False, ... data lies strictly below the bin
        10% percentile is calculated as: value at least 10% data is LOWER than
        """
        axis = self.get_axis_number(axis)
        
        # Shape of histogram
        s = self.histogram.shape

        # Shape of histogram after axis has been collapsed to 1
        s_collapsed = list(s)
        s_collapsed[axis] = 1

        # Shape of histogram with axis removed entirely
        s_removed = np.concatenate([s[:axis], s[axis + 1:]]).astype(np.int)

        # Using np.where here is too tricky, as it may not return a value for each "bin-columns"
        # First, get an array which has a minimum at the percentile-containing bins
        # The minimum may not be unique: if later bins are empty, they will not be 
        if inclusive:
            ecdf = self.cumulative_density(axis).histogram
        else:
            density = self.normalize(axis).histogram
            ecdf = ecdf - density
        ecdf = np.nan_to_num(ecdf)    # Since we're relying on self-equality later
        x = ecdf - 2 * (ecdf >= percentile / 100)
                
        # We now want to get the location of the minimum
        # To ensure it is unique, add a very very very small monotonously increasing bit to x
        # Nobody will want 1e-9th percentiles, right? TODO
        sz = np.ones(len(s), dtype=np.int)
        sz[axis] = -1
        x += np.linspace(0, 1e-9, s[axis]).reshape(sz)
        
        # 1. Find the minimum along the axis
        # 2. Reshape to s_collapsed and perform == to get a mask
        # 3. Apply the mask to the bin centers along axis
        # 4. Unflatten with reshape
        result = self.all_axis_bin_centers(axis)[
            x == np.min(x, axis=axis).reshape(s_collapsed)
        ]
        result = result.reshape(s_removed)

        if self.dimensions == 2:
            new_hist = Hist1d
        else:
            new_hist = Histdd
        return new_hist.from_histogram(histogram=result,
                                       bin_edges=itemgetter(*self.other_axes(axis))(self.bin_edges),
                                       axis_names=self.axis_names_without(axis))