如何在 Python 中计算给定百分位数内的平均值？

Question

我正在做一些科学计算，但我找不到执行以下操作的优雅方式。假设我有一个二维numpy 数组D，它在一天中多次存储给定数量的测量值。每行对应不同的测量仪器，每列对应一天中完成测量的不同时刻。

考虑所需百分位数的列表。例如：

quantiles = [0.25, 0.5, 0.75]

我的目标是按百分位组计算一天中每个时刻的平均测量值。换句话说，给定一列测量值，我想将该列中的所有测量值按照上述分位数分组，然后在组内取平均值。使用该示例，我将在一天中的每个时刻有 4 个组：下四分位数的测量值，然后是第 25 个和第 50 个四分位数之间的测量值，第 50 个和第 75 个四分位数之间的测量值，最后是最后一个四分位数的测量值。因此，如果m 是一天中进行测量的时刻数，q 是quantiles 变量中的元素数，我想要的输出将是qxm numpy 数组。

目前，我正在以最低效和最硬编码的方式执行此操作。我们开始：

quantiles = [0.25, 0.5, 0.75]
window = "30min"
moments = pd.date_range(start = "9:30", end = "16:00", freq = window).time
quantile_curves = np.zeros((len(quantiles)+1, len(moments)-1))
EmpQuantiles = np.quantile(D, quantiles, axis = 0)
for moment in range(len(moments)-1):
    quantile_curves[0, moment] = np.mean(D[:, moment][D[:,moment] < EmpQuantiles[0, moment]])
    quantile_curves[1, moment] = np.mean(D[:, moment][np.logical_and(D[:,moment] > EmpQuantiles[0, moment], D[:,moment] <EmpQuantiles[1, moment])])
    quantile_curves[2, moment] = np.mean(D[:, moment][np.logical_and(D[:,moment] > EmpQuantiles[1, moment], D[:,moment] <EmpQuantiles[2, moment])])
    quantile_curves[3, moment] = np.mean(D[:, moment][D[:,moment] > EmpQuantiles[2, moment]])

有什么优雅和简单的方法来做到这一点？我在这里找不到答案，但是R 中有一个相关（但不相同）的问题：ddply multiple quantiles by group

我打算绘制一天内组内平均值的演变。我展示了我在下面得到的情节（我对情节感到满意，并且得到了我想要的结果，但是我寻求更好的方法来计算 quantile_curves 变量）：

提前非常感谢！

Answer 1

您可以使用masked_arrays 有效地做到这一点：

import numpy as np

quantiles = [0.25, 0.5, 0.75]
print('quantiles:\n', quantiles)

moments = [f'moment {i}' for i in range(5)]
print('nb of moments:\n', len(moments))
nb_measurements = 10000
D = np.random.rand(nb_measurements,len(moments))
quantile_values = np.quantile(D,quantiles,axis=0)
print('quantile_values (for each moment):\n', quantile_values)

quantile_curves = np.zeros((len(quantiles)+1,len(moments)))
quantile_curves[0, :] = np.mean(np.ma.masked_array(D, mask=D>quantile_values[[0],:]), axis=0)
for q in range(len(quantiles)-1):
  quantile_curves[q+1, :] = np.mean(np.ma.masked_array(D, mask=np.logical_or(D<quantile_values[[q],:], D>quantile_values[[q+1],:])), axis=0)
quantile_curves[len(quantiles), :] = np.mean(np.ma.masked_array(D, mask=D<quantile_values[[len(quantiles)-1],:]), axis=0)

print('mean for each group and at each moment:')
print(quantile_curves)

输出：

% python3 script.py
quantiles:
 [0.25, 0.5, 0.75]
nb of moments:
 5
quantile_values (for each moment):
 [[0.25271343 0.25434056 0.24658732 0.24612319 0.25221014]
 [0.51114344 0.50103699 0.49671249 0.49113293 0.49819521]
 [0.75629377 0.75427293 0.74676209 0.74211813 0.7490436 ]]
mean for each group and at each moment
[[0.12650993 0.12823392 0.12492136 0.12200609 0.12655318]
 [0.3826476  0.373516   0.37050513 0.36974876 0.37722219]
 [0.63454102 0.63023986 0.62280545 0.61696283 0.6238492 ]
 [0.87866019 0.87614489 0.87492553 0.87253142 0.87403426]]

请注意，我使用的是 0 到 1 之间的随机值，这就是为什么分位数值（组区间的末端）几乎等于分位数的原因。也不是此代码适用于任意数量的分位数或矩。