使用NumPy或SciPy计算加权中位数。

Question

使用NumPy或SciPy计算加权中位数。

pythonnumpyscipystatisticssas-jmp

13

我正在尝试自动化 JMP 执行的一个过程（分析-> 分布，将列A输入为“Y值”，使用后续列作为“权重”值）。在JMP中，您必须一次处理一列 - 我想使用Python循环遍历所有列，并创建一个数组，显示每个列的中位数。

例如，如果质量数组为[0, 10, 20, 30]，并且第1列的权重数组为[30、191、9、0]，则质量数组的加权中位数应为10。然而，我不确定如何得出这个答案。

到目前为止，我已经:

导入了csv文件，将权重显示为数组，同时掩盖了值为0的数据；并
创建了一个与权重数组相同大小和形状的“Y值”数组(113x32)。虽然我不完全确定需要这样做，但认为它比使用for循环更容易进行加权。

我不确定接下来该怎么做。基本上，“Y值”是一系列质量值，数组中的所有列表示找到每个质量的数据点的数量。我需要根据报告频率找到中位数质量。

我不是 Python 或统计方面的专家，所以如果我省略了任何有用的细节，请让我知道！

更新：以下是我目前所做的一些代码：

#Boilerplate & Import files
import csv
import scipy as sp
from scipy import stats
from scipy.stats import norm
import numpy as np
from numpy import genfromtxt
import pandas as pd
import matplotlib.pyplot as plt

inputFile = '/Users/cl/prov.csv'
origArray = genfromtxt(inputFile, delimiter = ",")
nArray = np.array(origArray)
dimensions = nArray.shape
shape = np.asarray(dimensions)

#Mask values ==0
maTest = np.ma.masked_equal(nArray,0)

#Create array of masses the same shape as the weights (nArray)
fieldLength = shape[0]
rowLength = shape[1]

for i in range (rowLength):
    createArr = np.arange(0, fieldLength*10, 10)
    nCreateArr = np.array(createArr)
    massArr.append(nCreateArr)
    nCreateArr = np.array(massArr)
nmassArr = nCreateArr.transpose()

- Car

一些示例输入/输出数据会很有帮助，同时尽量展示你已经完成的代码。 - M4rtini

5个回答

8

如果我正确理解了您的问题，我们可以做的是总结观察结果，将其除以2得到对应于中位数的观察数量。从那里开始，我们需要弄清楚这个数字是哪个观察结果。

在这里的一个技巧是使用np.cumsum计算观察总和。它给出了一个运行累积总和。

例如： np.cumsum([1,2,3,4]) -> [1, 3, 6, 10] 每个元素都是所有先前元素和它本身的总和。我们有10个观察结果。因此平均值将是第5个观察结果。 (我们通过将最后一个元素除以2来得到5)。现在看一下cumsum的结果，我们可以很容易地看出这必须是第二个和第三个元素之间的观察结果(第3和第6个观察结果)。

因此，我们需要做的就是找出中位数(5)适合的索引位置。 np.searchsorted正好符合我们的需求。它会找到要插入一个元素以使数组保持排序的索引。

像这样编写代码：

import numpy as np
#my test data
freq_count = np.array([[30, 191, 9, 0], [10, 20, 300, 10], [10,20,30,40], [100,10,10,10], [1,1,1,100]])

c = np.cumsum(freq_count, axis=1) 
indices = [np.searchsorted(row, row[-1]/2.0) for row in c]
masses = [i * 10 for i in indices] #Correct if the masses are indeed 0, 10, 20,...

#This is just for explanation.
print "median masses is:",  masses
print freq_count
print np.hstack((c, c[:, -1, np.newaxis]/2.0))

输出将是：

median masses is: [10 20 20  0 30]  
[[ 30 191   9   0]  <- The test data
 [ 10  20 300  10]  
 [ 10  20  30  40]  
 [100  10  10  10]  
 [  1   1   1 100]]  
[[  30.   221.   230.   230.   115. ]  <- cumsum results with median added to the end.
 [  10.    30.   330.   340.   170. ]     you can see from this where they fit in.
 [  10.    30.    60.   100.    50. ]  
 [ 100.   110.   120.   130.    65. ]  
 [   1.     2.     3.   103.    51.5]]

- M4rtini

非常感谢您的解释！我已经接近了，但还没有完全理解。我认为我没有完全表达清楚我的问题——基本上，中位数应该始终是质量范围内的数字——[30、191、9、0]的频率分别对应于质量[0、10、20、30]（即质量在0-10范围内出现了30次，质量在10-20范围内出现了191次，等等）。根据您上面的答案，看起来我得到的是频率计数的中位数，对吗？ - Car

是的，它找到了频率计数的中位数，然后将其与质量相关联。利用这一点，质量范围直接与频率计数的元素相关。您需要找出真正的中位数还是包含中位数的范围？这将找到包含中位数的范围。 - M4rtini

你能否尝试提供更多的输入和输出示例，或者检查我使用的“测试数据”，并说明它们的输出应该是什么。 - M4rtini

理想情况下，我会找到真正的中位数，但范围也可以。使用您的测试数据，我分别找到了[20、25、25、25、25]的中位数。这是一些实际数据[30、191、30、0、0、0、0、0、0、0、0、0、0、0]，[0、99、256、254、82、5、0、0、0、0、0、0、0、0]，[0、0、0、65、205、189、249、120、72、40、2、0、0、0]，[0、0、0、0、0、1、59、192、324、204、188、127、104、29]。这些与从0-130计数的质量相对应，每10个为一组。使用JMP的中位数：[10、30、65、90]。 - Car

你编辑后的中位数是[125.5、348、471和614]。看起来已经接近了——它们逐渐变大，这与JMP的模式相同。我会继续调整它，看看是否有一个小的调整可以让它完成剩下的部分，但如果你有更多的建议，我会很感激！乍一看，可能是指数公式的问题——我得到的输出是0、10、50、80（在修改为(i-1)*10以从0开始时）。 - Car

6

wquantiles 是一个小型的 Python 包，它可以准确地完成你所需的功能。它在内部使用 np.cumsum() 和 np.interp()。

- muzzle

我已经有了加权 introselect 的适当实现计划好几年了 :( - Mad Physicist

1

我最终根据@muzzle和@maesers的回复编写了该函数:

def weighted_quantiles(values, weights, quantiles=0.5, interpolate=False):

    i = values.argsort()
    sorted_weights = weights[i]
    sorted_values = values[i]
    Sn = sorted_weights.cumsum()

    if interpolate:
        Pn = (Sn - sorted_weights/2 ) / Sn[-1]
        return np.interp(quantiles, Pn, sorted_values)
    else:
        return sorted_values[np.searchsorted(Sn, quantiles * Sn[-1])]

“interpolate True”和“interpolate False”的区别如下：

weighted_quantiles(np.array([1, 2, 3, 4]), np.ones(4))
> 2 
weighted_quantiles(np.array([1, 2, 3, 4]), np.ones(4), interpolate=True)
> 2.5

（对于如[1,2,3,4,5]这样的不均匀数组，没有区别）

速度测试显示，在未插值的情况下，它与@maesers的函数具有相同的性能，并且在插值的情况下它的性能是该函数的两倍。

- Mahé

0

分享一些我得到帮助的代码。这可以让您对Excel电子表格的每一列运行统计数据。

import xlrd
import sys
import csv
import numpy as np
import itertools
from itertools import chain

book = xlrd.open_workbook('/filepath/workbook.xlsx')
sh = book.sheet_by_name("Sheet1")
ofile = '/outputfilepath/workbook.csv'

masses = sh.col_values(0, start_rowx=1)  # first column has mass
age = sh.row_values(0, start_colx=1)   # first row has age ranges

count = 1
mass = []
for a in ages:
    age.append(sh.col_values(count, start_rowx=1))
    count += 1

stats = []
count = 0    
for a in ages:
    expanded = []
    # create a tuple with the mass vector

    age_mass = zip(masses, age[count])
    count += 1
    # replicate element[0] for element[1] times
    expanded = list(list(itertools.repeat(am[0], int(am[1]))) for am in age_mass)

    #  separate into one big list
    medianlist = [x for t in expanded for x in t]

    # convert to array and mask out zeroes
    npa = np.array(medianlist)
    npa = np.ma.masked_equal(npa,0)

    median = np.median(npa)
    meanMass = np.average(npa)
    maxMass = np.max(npa)
    minMass = np.min(npa)
    stdev = np.std(npa)

    stats1 = [median, meanMass, maxMass, minMass, stdev]
    print stats1

    stats.append(stats1)

np.savetxt(ofile, (stats), fmt="%d")

- Car

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

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

- masaers · Accepted Answer

既然这是在NumPy中加权中位数的谷歌搜索结果中排名最高的，我将添加一个简单的函数，以选择两个数组的加权中位数，而不改变它们的内容，并且对值的顺序没有任何假设（以防其他人也在此寻找相同的前提条件的快速方法）。

def weighted_median(values, weights):
    i = np.argsort(values)
    c = np.cumsum(weights[i])
    return values[i[np.searchsorted(c, 0.5 * c[-1])]]

使用argsort可以在不更改或复制它们的内容的情况下保持两个数组之间的对齐。很容易将其扩展到任意数量的任意分位数。

更新

由于很可能一开始不太明显如何轻松地扩展到任意分位数，因此这里有代码:

def weighted_quantiles(values, weights, quantiles=0.5):
    i = np.argsort(values)
    c = np.cumsum(weights[i])
    return values[i[np.searchsorted(c, np.array(quantiles) * c[-1])]]

默认为中位数，但可以传递任何分位数或分位数列表。返回类型与传入的quantiles相同，其中列表升级为NumPy数组。足够均匀分布的值，确实可以近似地表示输入：

>>> weighted_quantiles(np.random.rand(10000), np.random.rand(10000), [0.01, 0.05, 0.25, 0.50, 0.75, 0.95, 0.99])
array([0.01235101, 0.05341077, 0.25355715, 0.50678338, 0.75697424,0.94962936, 0.98980785])
>>> weighted_quantiles(np.random.rand(10000), np.random.rand(10000), 0.5)
0.5036283072043176
>>> weighted_quantiles(np.random.rand(10000), np.random.rand(10000), [0.5])
array([0.49851076])

更新2

在小数据集中，如果未实际观察到中位数/分位数，则能够在两个观测值之间插值可能非常重要。如果权重质量在它们之间平均分配（或者在分位点/1-分位点处），则可以通过计算两个数字之间的中点来很容易地添加此功能。由于需要条件，因此即使quantiles是单个标量，该函数也始终返回NumPy数组。现在输入也需要是NumPy数组（除了quantiles仍然可以是单个数字）。

def weighted_quantiles_interpolate(values, weights, quantiles=0.5):
    i = np.argsort(values)
    c = np.cumsum(weights[i])
    q = np.searchsorted(c, quantiles * c[-1])
    return np.where(c[q]/c[-1] == quantiles, 0.5 * (values[i[q]] + values[i[q+1]]), values[i[q]])

这个函数对于小于2的数组会失败（原始函数可以处理非空数组）。

>>> weighted_quantiles_interpolate(np.array([2, 1]), np.array([1, 1]), 0.5)
array(1.5)

请注意，在实际数据集中工作时，很少需要使用此扩展，因为我们通常具有 (a) 大型数据集和 (b) 实值权重，这使得最终准确地落在分位点边缘的可能性非常小，如果确实发生了这种情况，那可能是由于四舍五入误差导致的。然而出于完备性考虑，仍然包含此扩展。