描述性统计中分位数的计算与实现

在数据分析任务中经常需要计算分位数，如50、75、90、95、99分位数。一般项目中，基于Python语言时可以使用numpy library的percentile()方法，基于Java语言时可以使用Apache的commons-math3包中的StatUtils.percentile()方法分别进行计算。但实际上，这两种实现的计算逻辑却有所不同。甚至，关于分位数的具体计算逻辑在网上检索到的表述也存在较大差异。作为分位数算法的经典实现，numpy及commons-math3的官方文档虽然都提供了计算逻辑的文字描述，但并不十分清晰，源码中又包含了许多性能优化的成分，也不能直观体现分位数的计算逻辑。

对此，本文抽取了两个开源库中分位数计算逻辑的核心代码，从而不难推理出分位数的计算公式。

首先是numpy中的实现方式。

def percentile(values: list, p: int) -> float:
    if not (0 < p <= 100):
        raise Exception("P 的取值区间为 (0, 100]")

    size = len(values)
    if size == 0:
        return nan

    asc_values = sorted(values)
    if size == 1 or p < 1:
        return asc_values[0]
    if p == 100:
        return asc_values[-1]

    index = (p / 100) * (size - 1)
    index_below = floor(index)
    index_above = floor(index) + 1

    # asc_values[index_below] + (index - floor(index)) * (asc_values[index_above] - asc_values[index_below])
    return (asc_values[index_below] * (1 - index + index_below)) + (asc_values[index_above] * (index - index_below))

然后是commons-math3中的实现方式。

public static double percentile(double[] values, double p) throws Exception {
    if (p <= 0 || p > 100) {
        throw new Exception("P 的取值区间为 (0, 100]");
    }

    int size = values.length;
    if (size == 0) {
        return Double.NaN;
    }

    double[] ascValues = Arrays.stream(values).sorted().toArray();
    if (size == 1 || p < 1) {
        return ascValues[0];
    }
    if (p == 100) {
        return ascValues[size - 1];
    }

    double index = (p / 100) * (size + 1);
    int index_below = max((int) floor(index) - 1, 0);
    int index_above = min((int) floor(index), size - 1);

    return ascValues[index_below] + (index - floor(index)) * (ascValues[index_above] - ascValues[index_below]);
}

可见，二者的区别仅在于分位点两侧数组下标的选取上，这也是分位数计算的核心。在计算结果上，commons-math3中的实现也可能会比numpy中的实现偏大一点。

大数据背景下，上述分位数的计算方法需要完成对超大集合的排序，既会消耗大量计算资源，又无法满足实时计算的需求。当原始集合中的数值存在大量重复时，可以将排序过程简化为对去重后的子集进行排序。如此，只需要记录子集中各元素的下标范围，便可以快速检索出对原始集合排序后的某个下标所对应的具体元素，从而达到计算分位数的目的。参考numpy中算法，基于Python语言的具体实现如下。

def values_group(values: list) -> (int, list):
    bucket_map = dict()
    for value in values:
        bucket_map[value] = bucket_map.get(value, 0) + 1

    accumulation = 0
    bucket_list = list()
    for value in sorted(bucket_map.keys()):
        accumulation += bucket_map.get(value)
        bucket_list.append((value, accumulation))

    return accumulation, bucket_list


def values_index(values: list, index: int) -> float:
    for value in values:
        if index < value[1]:
            return value[0]


def percentile(values: list, p: int) -> float:
    if not (0 < p <= 100):
        raise Exception("P 的取值区间为 (0, 100]")

    size, asc_values = values_group(values)
    if size == 0:
        return nan

    if size == 1 or p < 1:
        return asc_values[0][0]
    if p == 100:
        return asc_values[-1][0]

    index = (p / 100) * (size - 1)
    index_below = floor(index)
    index_above = floor(index) + 1

    return (values_index(asc_values, index_below) * (1 - index + index_below)) + (
            values_index(asc_values, index_above) * (index - index_below))

上述values_group()方法可以借用Spark等大数据计算引擎的groupBy及开窗算子进行计算。