Faster mquantiles function - simplified

gistfile1.txt

def compute_quantile_curve(curves, quantile):
    """
    Compute the quantile aggregate of a set of curves. This method is used in
    the case where hazard curves are computed using the Monte-Carlo logic tree
    sampling approach. In this case, the weights are implicit.

    :param curves:
        2D array-like collection of hazard curve PoE values. Each element
        should be a sequence of PoE `float` values. Example::

            [[0.5, 0.4, 0.3], [0.6, 0.59, 0.1]]
    :param float quantile:
        The quantile value. We expected a value in the range [0.0, 1.0].

    :returns:
        A numpy array representing the quantile aggregate of the input
        ``curves`` and ``quantile``.
    """
    # this implementation is an alternative to:
    # return numpy.array(mstats.mquantiles(curves, prob=quantile, axis=0))[0]

    # more or less copied from the scipy mquantiles function, just faster

    arr = numpy.array(curves)

    p = numpy.array(quantile)
    m = 0.4 + p * 0.2

    n = len(arr)
    aleph = n * p + m
    k = numpy.floor(aleph.clip(1, n - 1)).astype(int)
    gamma = (aleph-k).clip(0, 1)

    data = numpy.sort(arr, axis=0).transpose()
    return (1.0 - gamma) * data[:, k-1] + gamma * data[:, k]