10

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

我正在尝试计算每个 1 分钟间隔值的平均值/标准差（取 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 分钟只是一个例子）

4

6 回答 6

11

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

行。我觉得有趣的是，对我来说只有一个合乎逻辑的解决方案，但许多其他人找到了其他解决方案。无论如何，解决方案似乎很简单。给定向量 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);
于 2010-02-24T11:17:06.670 回答
4

您可以尝试创建一个元胞数组并通过 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]），并且由于间隔的开始不完全相同.

于 2010-02-24T02:25:31.797 回答
3

这是一种使用二分搜索的方法。9900 个元素的速度提高了 6-10 倍，99900 个元素的速度提高了大约 64 倍。仅使用 900 个元素很难获得可靠的时间，所以我不确定在那个尺寸下哪个更快。如果您考虑直接从生成的数据进行 tx，它几乎不会使用额外的内存。除此之外，它只有四个额外的浮点变量（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) 以你正在做的方式找到它们。

于 2010-02-24T05:40:30.800 回答
2

您可以indices使用 bsxfun 一次计算所有内容：

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

这比循环更快，但需要一次存储它们（时间与空间的权衡）..

于 2010-02-24T04:11:59.750 回答
2

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

您可以通过自己进行一些棘手的累积和、索引以及计算均值和标准差来避免循环或使用元胞数组。这里有一些我相信会起作用的代码，尽管我不确定它如何在速度方面与其他解决方案相提并论：

[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

上面使用此 Wikipedia 页面上的公式的简化计算标准偏差。

于 2010-02-24T06:31:43.870 回答
0

与上述相同的答案，但具有参数区间 ( 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);
于 2013-12-02T14:37:28.120 回答