Python:使用梯形规则快速计算平均值

Question

使用Python，我必须处理一些数据。

我得到了大约 50 个函数在大约 1000 万个时间点的值。这些值以二维列表矩阵的形式给出，即matrix[i]是一个值列表

[t_i，t_i 处 f1 的值，t_i 处 f2 的值 ...，t_i 处 fN 的值]

其中N = 50。

由于

• 数据中可能存在的噪声(函数是一些测量值)
• 非等距时间点(有时时间步长是几秒，但有时可以以天为单位)

我决定使用给定值在固定长度的某些预定义时间间隔上的平均值。

我尝试了不同长度的间隔:范围在一分钟到一小时之间。

我计算平均值的算法如下:

matrix = ...  # read matrix
t0 = matrix[0][0]
ts_new = [t0 + i * time_step for i in range(some_bound)]
buckets = [[] for t in ts_new]
for i, instance in enumerate(matrix):
    t_i = instance[0]
    put i to the bucket with index j, such that ts_new[j] <= t_i < ts_new[j + 1] 
for bucket in buckets:
    compute the average values of f1, ... , fN over the instances from bucket

The bottle-neck of the algorithm is the last for-loop.

如果我将 matrix 转换为 numpy.array 并将 bucket 的平均值计算为 matrix[bucket, :] .mean(axis=0)，这工作得相当快，但计算值没有多大意义:

如果f1在ts = [0, 99, 100]时刻的值分别为ys = [0, 0, 2]，mean 函数返回 2/3(如预期)。但是，f1 的平均值应该更接近于0。使用梯形规则，可以得到平均值0.01，这样更有意义。

所以，目前，我正在使用

• 计算桶平均值的scipy.integrate.trapz方法:梯形面积除以间隔长度
• scipy.interpolate.interp1d 方法用于获取函数 f 在间隔边界处的值，例如，我使用前一个桶中的最后一个点和第一个点给定存储桶中的点，计算相应时间间隔开始时的值(时间间隔结束时的值类似)

需要第二个项目符号，因为一分钟的时间间隔长度非常短，有时桶中只有一两个点。程序是这样的:

# for one bucket
means = [0 for col in range(N)]
for col in range(1, N + 1):  # for each function f
    xs = []
    ys = []
    if can_interpolate_at_start:
        f_lin = scipy.interpolate.interp1d([tPrevLast, tNowFirst], [yPrevLast, yNowFirst])
        xs.append(t_bucketStart)
        ys.append(f_lin(t_bucketStart))
    xs += matrix[bucket, 0]
    ys += matrix[bucket, col]
    if can_interpolate_at_end:
        # ...
    means[col - 1] = scipy.integrate.trapz(ys, xs) / (xs[-1] - xs[0])

can_interpolate_at_start 和 can_interpolate_at_end 的值仅取决于时间间隙(但必须特别注意第一个和最后一个存储桶...):我不如果前一个存储桶中的最后一个点与当前存储桶中的时间差太大，则使用插值点。

My problem: the current approach is really slow (two hours or so for one-minute time intervals). How can I make it faster?

Answer 1

您可以做的一个简单的事情是，从间隔的角度考虑，取每个间隔的每个函数的平均值，然后将 if 乘以间隔长度，然后除以总时间:

import numpy as np

matrix = ...
data = np.asarray(matrix)
t_diff = np.diff(data[:, 0])
means_sum = np.sum(t_diff[:, np.newaxis] * (data[:-1, 1:] + data[1:, 1:]) / 2, axis=0)
means = means_sum / (data[-1, 0] - data[0, 0])