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

Question

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

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

quantiles = [0.25, 0.5, 0.75]

我的目标是在一天中的每一时刻，按百分比组计算平均测量值。换句话说，给定一列度量值，我想根据上面的分位数对该列中的所有度量值进行分组排序，然后在分组中取平均值。使用这个例子，我会在一天中的每个时刻有4组：下四分位的测量值，然后是第25和第50个四分位之间的测量值，第50和第75个四分位之间的测量值，最后是最后一个四分位中的测量值。因此，如果m是进行测量的一天中的时刻数，而q是quantiles变量中的元素数，则我所需的输出将是qxmnumpy数组

目前，我正在以最低效和最硬编码的方式进行这项工作。我们开始：

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之间的随机值，这就是为什么分位数值（组间隔的极值）几乎等于分位数。也不是说这段代码适用于任意数量的分位数或矩