在 Numpy 中重新创建 R 分位数类型 2

Question

我正在将一些遗留代码从 R 迁移到 Python，但我无法匹配 quantile results with numpy percentile。

给定以下数字列表：

a1 = [
    5.75,6.13333333333333,7.13636363636364,9,10.1,4.80952380952381,8.82926829268293,4.7906976744186,3.83333333333333,6,6.1,
    8.88235294117647,30,5.7,3.98507462686567,6.83333333333333,8.39805825242718,4.78260869565217,7.26356589147287,5.67857142857143,
    3.58333333333333,6.69230769230769,14.3333333333333,14.3333333333333,5.125,5.16216216216216,5.36363636363636,10.7142857142857,
    4.90909090909091,7.5,8,6,6.93939393939394,10.4,6,6.8,5.33333333333333,10.3076923076923,4.5625,5.4,6.44,3.36363636363636,
    11.1666666666667,4.5,7.35714285714286,10.6363636363636,9.26746031746032,3.83333333333333,5.75,9.14285714285714,8.27272727272727,
    5,5.92307692307692,5.23076923076923,4.09375,6.25,4.63888888888889,6.07142857142857,5,5.42222222222222,3.93892045454545,4.8,
    8.71428571428571,6.25925925925926,4.12,5.30769230769231,4.26086956521739,5.22222222222222,4.64285714285714,5,3.64705882352941,
    5.33333333333333,3.65217391304348,3.54166666666667,10.0952380952381,3.38235294117647,8.67123287671233,2.66666666666667,3.5,4.875,
    4.5,6.2,5.45454545454545,4.89189189189189,4.71428571428571,1,5.33333333333333,6.09090909090909,4.36756756756757,6,5.17197452229299,
    4.48717948717949,5.01219512195122,4.83098591549296,5.25,8.52,5.47692307692308,5.45454545454545,8.6578947368421,8.35714285714286,3.25,
    8.5,4,5.95652173913043,7.05882352941176,7.5,8.6,8.49122807017544,5.14285714285714,4,13.3294117647059,9.55172413793103,5.57446808510638,
    4.5,8,4.11764705882353,3.9,5.14285714285714,6,4.66666666666667,6,3.75,4.93333333333333,4.5,5.21666666666667,6.53125,6,7,7.28333333333333,
    7.34615384615385,7.15277777777778,8.07936507936508,11.609756097561
]

在 R 中使用分位数使得

quantile(a1, probs=.05, type=2)

给出 3.541667

的结果

尝试 numpy 中的所有插值方法以找到相同的结果：

{x:np.percentile(a1,q=5, interpolation=x) for x in ['linear','lower','higher','nearest','midpoint']}

产量

{'linear': 3.566666666666666,
 'lower': 3.54166666666667,
 'higher': 3.58333333333333,
 'nearest': 3.58333333333333,
 'midpoint': 3.5625}

我们可以看到lower插值方法returns与R分位数类型2相同的结果

然而，在 R 中再次使用不同的分位数，我们得到不同的结果：

quantile(a1, probs=.95, type=2)

给出 10.71429

的结果

和 numpy:

{x:np.percentile(a1,q=95, interpolation=x) for x in ['linear','lower','higher','nearest','midpoint']}

产量

{'linear': 10.667532467532439,
 'lower': 10.6363636363636,
 'higher': 10.7142857142857,
 'nearest': 10.6363636363636,
 'midpoint': 10.67532467532465}

本例higher插值法returns结果相同

我希望足够熟悉 w/the R 分位数类型的人可以帮助我在 numpy 中重现相同的分位数逻辑。

Answer 1

您可以自己实现。使用 type=2 这是一个相当简单的计算。您要么采用下一个最高阶统计数据，要么采用不连续性（即 100 个值，并且您希望 p=0.06，正好落在第 6 个值上），您采用该阶统计数据和下一个最大阶统计数据的平均值。

import numpy as np

def R_type2(arr, p):
    """
    arr : array-like
    p : float between [0, 1]
    """
    #m=0 for Q_2(p) in R
    x = np.sort(arr)
    n = len(x)

    aleph = n*p
    k = np.floor(np.array(aleph).clip(1, n-1)).astype(int)
    gamma = {False: 1, True: 0.5}.get(aleph==k)   # Discontinuity or not

    # Deal with case where it should be smallest value 
    if aleph < 1:
        return x[k-1]  # x[0]
    else:
        return (1.-gamma)*x[k-1] + gamma*x[k]


R_type2(a1, 0.05)
#3.54166666666667

R_type2(a1, 0.95)
#10.7142857142857

提醒一句。 k 将是一个整数，而 n*p 是一个浮点数。通常，执行 aleph==k 是一个非常糟糕的主意，因为这会导致浮点数不准确的问题。。但是，因为 R 似乎实现了纯粹的相等性检查，为了保持一致性，我像上面那样保留了它。

我个人更喜欢从平等转变：{False: 1, True: 0.5}.get(aleph==k) {False: 1, True: 0.5}.get(np.isclose(aleph,k)) 这样浮点数问题就不会成为问题。

在 Numpy 中重新创建 R 分位数类型 2

Recreating R Quantile Type 2 in Numpy

python

r

pandas

numpy

quantile