使用numpy计算加权百分位数

Question

使用numpy计算加权百分位数

pythonnumpyweightedpercentile

53

有没有一种方法可以使用numpy.percentile函数计算加权百分位数？或者是否有其他的Python函数可用于计算加权百分位数？

谢谢！

- user308827

在我看来，Sam A 在下面的解决方案（https://dev59.com/72Eh5IYBdhLWcg3w9XWn#63440143）似乎是当前最佳实践的有力竞争者。 - geotheory

12个回答

15

这似乎现在已经在statsmodels中实现了。

from statsmodels.stats.weightstats import DescrStatsW
wq = DescrStatsW(data=np.array([1, 2, 9, 3.2, 4]), weights=np.array([0.0, 0.5, 1.0, 0.3, 0.5]))
wq.quantile(probs=np.array([0.1, 0.9]), return_pandas=False)
# array([2., 9.])

DescrStatsW对象还实现了其他方法，例如加权平均值等。 https://www.statsmodels.org/stable/generated/statsmodels.stats.weightstats.DescrStatsW.html

- Sam A.

14

一种快速的解决方法是先排序，然后进行插值：

def weighted_percentile(data, percents, weights=None):
    ''' percents in units of 1%
        weights specifies the frequency (count) of data.
    '''
    if weights is None:
        return np.percentile(data, percents)
    ind=np.argsort(data)
    d=data[ind]
    w=weights[ind]
    p=1.*w.cumsum()/w.sum()*100
    y=np.interp(percents, p, d)
    return y

- Kambrian

5

weighted_percentile(np.array([0,3,6,9]),50,weights=np.array([1,3,3,1])) 和 weighted_percentile(np.array([0,3,3,3,6,6,6,9]),50,weights=None) 的结果不同。 - Peter9192

将1.*w.cumsum()更改为(1.*w.cumsum()-0.5*w)（根据@imbr的回答），可以得到上述期望的结果。 - undefined

13

使用这个参考，可以更加简洁清晰地使用加权百分位数法。

import numpy as np

def weighted_percentile(data, weights, perc):
    """
    perc : percentile in [0-1]!
    """
    ix = np.argsort(data)
    data = data[ix] # sort data
    weights = weights[ix] # sort weights
    cdf = (np.cumsum(weights) - 0.5 * weights) / np.sum(weights) # 'like' a CDF function
    return np.interp(perc, cdf, data)

- imbr

11

我不知道加权百分位是什么意思，但从@Joan Smith的回答中看来，你只需要重复ar中的每个元素，可以使用numpy.repeat()函数：

我不熟悉“加权分位数”的含义，但根据 @Joan Smith 的回答，似乎您只需要在 ar 中重复每个元素，您可以使用 numpy.repeat() 函数：

import numpy as np
np.repeat([1,2,3], [4,5,6])

结果为：

array([1, 1, 1, 1, 2, 2, 2, 2, 2, 3, 3, 3, 3, 3, 3])

- HYRY

2

我认为这是更好的（更有效率的）答案。 - FooBar

18

然而，这仅支持整数权重。对于较大的数据集，很可能会占用大量内存和CPU时间。 - PiHalbe

3

从个人经验来看，我可以确认这种方法绝对不高效。如果您的向量很长且权重很大，您的计算机很快就会达到内存限制。 - geotheory

那也不适用于非整数权重。 - amyrit

10

对于额外的（非原创）答案表示歉意（没有足够的声誉在@nayyarv上发表评论）。他的解决方案对我有效（即它复制了np.percentage的默认行为），但我认为您可以通过查看原始np.percentage的编写方式来消除for循环。

def weighted_percentile(a, q=np.array([75, 25]), w=None):
    """
    Calculates percentiles associated with a (possibly weighted) array

    Parameters
    ----------
    a : array-like
        The input array from which to calculate percents
    q : array-like
        The percentiles to calculate (0.0 - 100.0)
    w : array-like, optional
        The weights to assign to values of a.  Equal weighting if None
        is specified

    Returns
    -------
    values : np.array
        The values associated with the specified percentiles.  
    """
    # Standardize and sort based on values in a
    q = np.array(q) / 100.0
    if w is None:
        w = np.ones(a.size)
    idx = np.argsort(a)
    a_sort = a[idx]
    w_sort = w[idx]

    # Get the cumulative sum of weights
    ecdf = np.cumsum(w_sort)

    # Find the percentile index positions associated with the percentiles
    p = q * (w.sum() - 1)

    # Find the bounding indices (both low and high)
    idx_low = np.searchsorted(ecdf, p, side='right')
    idx_high = np.searchsorted(ecdf, p + 1, side='right')
    idx_high[idx_high > ecdf.size - 1] = ecdf.size - 1

    # Calculate the weights 
    weights_high = p - np.floor(p)
    weights_low = 1.0 - weights_high

    # Extract the low/high indexes and multiply by the corresponding weights
    x1 = np.take(a_sort, idx_low) * weights_low
    x2 = np.take(a_sort, idx_high) * weights_high

    # Return the average
    return np.add(x1, x2)

# Sample data
a = np.array([1.0, 2.0, 9.0, 3.2, 4.0], dtype=np.float)
w = np.array([2.0, 1.0, 3.0, 4.0, 1.0], dtype=np.float)

# Make an unweighted "copy" of a for testing
a2 = np.repeat(a, w.astype(np.int))

# Tests with different percentiles chosen
q1 = np.linspace(0.0, 100.0, 11)
q2 = np.linspace(5.0, 95.0, 10)
q3 = np.linspace(4.0, 94.0, 10)
for q in (q1, q2, q3):
    assert np.all(weighted_percentile(a, q, w) == np.percentile(a2, q))

- grovduck

这很有用。但是，我不得不将 idx_high[idx_high > ecdf.size - 1] = ecdf.size - 1 包装在条件语句中，以使其适用于单个百分位数。猜想这就是为什么 numpy 源代码中有 zerod 的原因。 - Peter9192

4

weightedcalcs 包支持分位数：

import weightedcalcs as wc
import pandas as pd

df = pd.DataFrame({'v': [1, 2, 3], 'w': [3, 2, 1]})
calc = wc.Calculator('w')  # w designates weight

calc.quantile(df, 'v', 0.5)
# 1.5

- Max Ghenis

3

如评论中所述，对于浮点权重来说，简单地重复值是不可能的，对于非常大的数据集来说也不实用。有一个库可以在这里执行加权百分位数计算： http://kochanski.org/gpk/code/speechresearch/gmisclib/gmisclib.weighted_percentile-module.html 它对我很有帮助。

- Qwerty

2

我使用这个函数来满足我的需求:

def quantile_at_values(values, population, weights=None):
    values = numpy.atleast_1d(values).astype(float)
    population = numpy.atleast_1d(population).astype(float)
    # if no weights are given, use equal weights
    if weights is None:
        weights = numpy.ones(population.shape).astype(float)
        normal = float(len(weights))
    # else, check weights                  
    else:                                           
        weights = numpy.atleast_1d(weights).astype(float)
        assert len(weights) == len(population)
        assert (weights >= 0).all()
        normal = numpy.sum(weights)                    
        assert normal > 0.
    quantiles = numpy.array([numpy.sum(weights[population <= value]) for value in values]) / normal
    assert (quantiles >= 0).all() and (quantiles <= 1).all()
    return quantiles

我已经进行了向量化处理。
它有很多合理性检查。
它使用浮点数作为权重。
它可以不使用权重（→等权重）来工作。
它可以同时计算多个分位数。

如果要得到百分位数而非分位数，请将结果乘以100。

- PiHalbe

1

请注意，此函数返回给定值处的分位数，虽然与之相关，但并未回答问题，该问题是关于百分位数的（而且百分位数不等于分位数乘以100）。 - ILoveCoding

2

def weighted_percentile(a, percentile = np.array([75, 25]), weights=None):
    """
    O(nlgn) implementation for weighted_percentile.
    """
    percentile = np.array(percentile)/100.0
    if weights is None:
        weights = np.ones(len(a))
    a_indsort = np.argsort(a)
    a_sort = a[a_indsort]
    weights_sort = weights[a_indsort]
    ecdf = np.cumsum(weights_sort)

    percentile_index_positions = percentile * (weights.sum()-1)+1
    # need the 1 offset at the end due to ecdf not starting at 0
    locations = np.searchsorted(ecdf, percentile_index_positions)

    out_percentiles = np.zeros(len(percentile_index_positions))

    for i, empiricalLocation in enumerate(locations):
        # iterate across the requested percentiles 
        if ecdf[empiricalLocation-1] == np.floor(percentile_index_positions[i]):
            # i.e. is the percentile in between 2 separate values
            uppWeight = percentile_index_positions[i] - ecdf[empiricalLocation-1]
            lowWeight = 1 - uppWeight

            out_percentiles[i] = a_sort[empiricalLocation-1] * lowWeight + \
                                 a_sort[empiricalLocation] * uppWeight
        else:
            # i.e. the percentile is entirely in one bin
            out_percentiles[i] = a_sort[empiricalLocation]

    return out_percentiles

这是我的函数，它提供与...相同的行为。

np.percentile(np.repeat(a, weights), percentile)

使用更少的内存开销。np.percentile是一个O(n)的实现，因此对于小权重来说它可能更快。它已经解决了所有边缘情况——这是一个精确的解决方案。上面的插值答案假定是线性的，但在大多数情况下都是步进的，除非权重为1。

假设我们有数据[1,2,3]和权重[3, 11, 7]，我想要25%的百分位数。我的ecdf将会是[3, 10, 21]，我要找的是第5个值。插值将会看到[3,1]和[10,2]作为匹配并进行插值，尽管完全处于第二个bin中，其值为2，但仍会得到1.28的结果。

- nayyarv

这是完全错误的！ - Yahya

网页内容由stack overflow 提供, 点击上面的

可以查看英文原文，
原文链接

- Alleo · Accepted Answer

完全矢量化的numpy解决方案

这是我使用的代码。虽然不是最优解(我无法使用numpy编写最优解)，但仍比已接受的解决方案更快、更可靠。

def weighted_quantile(values, quantiles, sample_weight=None, 
                      values_sorted=False, old_style=False):
    """ Very close to numpy.percentile, but supports weights.
    NOTE: quantiles should be in [0, 1]!
    :param values: numpy.array with data
    :param quantiles: array-like with many quantiles needed
    :param sample_weight: array-like of the same length as `array`
    :param values_sorted: bool, if True, then will avoid sorting of
        initial array
    :param old_style: if True, will correct output to be consistent
        with numpy.percentile.
    :return: numpy.array with computed quantiles.
    """
    values = np.array(values)
    quantiles = np.array(quantiles)
    if sample_weight is None:
        sample_weight = np.ones(len(values))
    sample_weight = np.array(sample_weight)
    assert np.all(quantiles >= 0) and np.all(quantiles <= 1), \
        'quantiles should be in [0, 1]'

    if not values_sorted:
        sorter = np.argsort(values)
        values = values[sorter]
        sample_weight = sample_weight[sorter]

    weighted_quantiles = np.cumsum(sample_weight) - 0.5 * sample_weight
    if old_style:
        # To be convenient with numpy.percentile
        weighted_quantiles -= weighted_quantiles[0]
        weighted_quantiles /= weighted_quantiles[-1]
    else:
        weighted_quantiles /= np.sum(sample_weight)
    return np.interp(quantiles, weighted_quantiles, values)

示例:

weighted_quantile([1, 2, 9, 3.2, 4], [0.0, 0.5, 1.])

数组 [ 1. , 3.2, 9. ]

weighted_quantile([1, 2, 9, 3.2, 4], [0.0, 0.5, 1.], sample_weight=[2, 1, 2, 4, 1])

数组 [ 1. , 3.2, 9. ]