MATLAB：计算时间序列的每个 1 分钟间隔的平均值

Question

我有一堆时间序列，每个时间序列由两个组件描述，一个时间戳向量（以秒为单位）和一个测量值向量。时间向量是不均匀的（即以非规则间隔采样）

我试图计算每个 1 分钟间隔值的平均值/SD（采用 X 分钟间隔，计算其平均值，采用下一个间隔，.. .)。

我当前的实现使用循环。这是我到目前为止所拥有的示例：

t = (100:999)' + rand(900,1);       %' non-uniform time
x = 5*rand(900,1) + 10;             % x(i) is the value at time t(i)

interval = 1;         % 1-min interval
tt = ( floor(t(1)):interval*60:ceil(t(end)) )';  %' stopping points of each interval
N = length(tt)-1;

mu = zeros(N,1);
sd = zeros(N,1);

for i=1:N
    indices = ( tt(i) <= t & t < tt(i+1) ); % find t between tt(i) and tt(i+1)
    mu(i) = mean( x(indices) );
    sd(i) = std( x(indices) );
end

我想知道是否有更快的矢量化解决方案。这很重要，因为我有大量时间序列需要处理，每个时间序列比上面显示的示例要长得多。

欢迎任何帮助。

---

感谢大家的反馈。

我更正了 t 的生成方式，使其始终单调递增（排序），这并不是真正的问题。

另外，我可能没有明确说明这一点，但我的目的是为任何问题提供解决方案间隔长度（以分钟为单位）（1 分钟只是一个示例）

Answer 1

唯一合乎逻辑的解决方案似乎是......

好吧。我觉得很有趣的是，对我来说只有一种合乎逻辑的解决方案，但许多其他人找到了其他解决方案。无论如何，解决方案看起来确实很简单。给定向量 x 和 t，以及一组等距断点 tt，

t = sort((100:999)' + 3*rand(900,1));     % non-uniform time
x = 5*rand(900,1) + 10;             % x(i) is the value at time t(i)

tt = ( floor(t(1)):1*60:ceil(t(end)) )';

（请注意，我在上面对 t 进行了排序。）

我将用三行完全矢量化的代码来完成此操作。首先，如果中断是任意的并且间距可能不相等，我将使用 histc 来确定数据系列属于哪个间隔。鉴于它们是均匀的，只需执行以下操作：

int = 1 + floor((t - t(1))/60);

同样，如果不知道 t 的元素已排序，我会使用 min(t) 而不是 t(1)。完成此操作后，使用 Accumarray 将结果减少为平均值和标准差。

mu = accumarray(int,x,[],@mean);
sd = accumarray(int,x,[],@std);

Answer 2

您可以尝试创建一个元胞数组并通过 cellfun 应用平均值和标准差。对于 900 个条目，它比您的解决方案慢约 10%，但对于 90000 个条目，速度快约 10 倍。

[t,sortIdx]=sort(t); %# we only need to sort in case t is not monotonously increasing
x = x(sortIdx);

tIdx = floor(t/60); %# convert seconds to minutes - can also convert to 5 mins by dividing by 300
tIdx = tIdx - min(tIdx) + 1; %# tIdx now is a vector of indices - i.e. it starts at 1, and should go like your iteration variable.

%# the next few commands are to count how many 1's 2's 3's etc are in tIdx
dt = [tIdx(2:end)-tIdx(1:end-1);1]; 
stepIdx = [0;find(dt>0)];
nIdx = stepIdx(2:end) - stepIdx(1:end-1); %# number of times each index appears

%# convert to cell array
xCell = mat2cell(x,nIdx,1);

%# use cellfun to calculate the mean and sd
mu(tIdx(stepIdx+1)) = cellfun(@mean,xCell); %# the indexing is like that since there may be missing steps
sd(tIdx(stepIdx+1)) = cellfun(@mean,xCell);

注意：我的解决方案不会给出与您的完全相同的结果，因为您在最后跳过了一些时间值（1:60:90 是 [1,61]），并且由于间隔的开始并不完全相同。

Answer 3

这是一种使用二分搜索的方法。对于 9900 个元素，速度提高了 6-10 倍；对于 99900 个元素，速度提高了约 64 倍。仅使用 900 个元素很难获得可靠的时间，因此我不确定在该尺寸下哪个更快。如果您考虑直接从生成的数据进行交易，它几乎不使用额外的内存。除此之外，它只有四个额外的浮点变量（prevind、first、mid 和 last）。

% Sort the data so that we can use binary search (takes O(N logN) time complexity).
tx = sortrows([t x]);

prevind = 1;

for i=1:N
    % First do a binary search to find the end of this section
    first = prevind;
    last = length(tx);
    while first ~= last
        mid = floor((first+last)/2);
        if tt(i+1) > tx(mid,1)
            first = mid+1;
        else
            last = mid;
        end;
    end;
    mu(i) = mean( tx(prevind:last-1,2) );
    sd(i) = std( tx(prevind:last-1,2) );
    prevind = last;
end;

它使用您最初拥有的所有变量。我希望它适合您的需求。它更快，因为使用二分搜索查找索引需要 O(log N) 时间，但按照您的方式查找索引需要 O(N) 时间。

Answer 4

您可以使用 bsxfun 一次性计算所有索引：

indices = ( bsxfun(@ge, t, tt(1:end-1)') & bsxfun(@lt, t, tt(2:end)') );

这比循环更快，但需要一次存储所有索引（时间与空间权衡）。

Answer 5

免责声明：我在纸上解决了这个问题，但还没有机会“在计算机中”检查它...

您可以通过执行一些棘手的累积来避免循环或使用元胞数组自己求和、索引并计算平均值和标准差。下面是一些我相信可以工作的代码，尽管我不确定它如何与其他解决方案相比提高速度：

[t,sortIndex] = sort(t);  %# Sort the time points
x = x(sortIndex);         %# Sort the data values
interval = 60;            %# Interval size, in seconds

intervalIndex = floor((t-t(1))./interval)+1;  %# Collect t into intervals
nIntervals = max(intervalIndex);              %# The number of intervals
mu = zeros(nIntervals,1);                     %# Preallocate mu
sd = zeros(nIntervals,1);                     %# Preallocate sd

sumIndex = [find(diff(intervalIndex)) ...
            numel(intervalIndex)];  %# Find indices of the interval ends
n = diff([0 sumIndex]);             %# Number of samples per interval
xSum = cumsum(x);                   %# Cumulative sum of x
xSum = diff([0 xSum(sumIndex)]);    %# Sum per interval
xxSum = cumsum(x.^2);               %# Cumulative sum of x^2
xxSum = diff([0 xxSum(sumIndex)]);  %# Squared sum per interval

intervalIndex = intervalIndex(sumIndex);  %# Find index into mu and sd
mu(intervalIndex) = xSum./n;                             %# Compute mean
sd(intervalIndex) = sqrt((xxSum-xSum.*xSum./n)./(n-1));  %# Compute std dev

上面使用 此维基百科页面上找到的公式的简化。

Answer 6

与上面相同的答案，但具有参数间隔（window_size）。
矢量长度的问题也得到解决。

window_size = 60; % but it can be any value 60 5 0.1, which wasn't described above

t = sort((100:999)' + 3*rand(900,1));     % non-uniform time
x = 5*rand(900,1) + 10;                   % x(i) is the value at time t(i)

int = 1 + floor((t - t(1))/window_size);
tt = ( floor(t(1)):window_size:ceil(t(end)) )';



% mean val and std dev of the accelerations at speed
mu = accumarray(int,x,[],@mean);
sd = accumarray(int,x,[],@std);

%resolving some issue with sizes (for i.e. window_size = 1 in stead of 60)
while ( sum(size(tt) > size(mu)) > 0 ) 
  tt(end)=[]; 
end

errorbar(tt,mu,sd);