我正在做一些科学计算,但我找不到一种优雅的方法来执行下面的操作。假设我有一个二维numpy
数组D
,它存储一天中几次给定量的测量值。每行对应于不同的测量仪器,每列对应于一天中完成测量的不同时刻。
考虑所需百分位数的列表。例如:
quantiles = [0.25, 0.5, 0.75]
我的目标是在一天中的每个时刻按百分位数组计算平均测量值。换句话说,给定一列测量值,我希望根据上面的分位数对该列中的所有测量值进行分组排序,然后在组内取平均值。使用这个例子,我将在一天中的每个时刻有4组:低四分位数的测量值,然后是第25和50个四分位数之间的测量值,50和75之间的测量值,最后是最后四分位数的测量值。因此,如果m
是一天中进行测量的时刻数,q
是quantiles
变量中的元素数,那么我想要的输出将是q
xm
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
变量):
提前谢谢你!
发布于 2020-07-11 03:09:17
您可以使用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之间的随机值,这就是为什么分位数(组间隔的末端)几乎等于分位数。同样,这段代码也不适用于任意数量的分位数或矩。
https://stackoverflow.com/questions/62839196
复制相似问题