MATLAB: вычислить среднее значение каждого 1-минутного интервала временного ряда

Question

У меня есть набор временных рядов, каждый из которых описывается двумя компонентами: вектором метки времени (в секундах) и вектором измеренных значений. Вектор времени является неоднородным (то есть дискретизируется с нерегулярными интервалами)

Я пытаюсь вычислить среднее значение / SD для каждого 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 минута была только примером)

Answer 1

Кажется, единственное логическое решение ...

Ok. Мне смешно, что для меня есть только одно логическое решение, но многие другие находят другие решения. В любом случае, решение кажется простым. При заданных векторах 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)) )';

(Обратите внимание, что я отсортировал т выше.)

Я бы сделал это в трех полностью векторизованных строках кода. Во-первых, если бы разрывы были произвольными и потенциально неравными в интервале, я бы использовал historyc, чтобы определить, в какие интервалы попадает ряд данных. Учитывая, что они однородны, просто сделайте это:

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. Это примерно на 10% медленнее, чем ваше решение для 900 записей, но примерно в 10 раз быстрее для 90000 записей.

[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

Вот способ, который использует бинарный поиск . Это в 6-10 раз быстрее для 9900 элементов и примерно в 64 раза быстрее для 99900 элементов. Трудно было получить надежное время, используя только 900 элементов, поэтому я не уверен, что быстрее при таком размере. Он почти не использует дополнительную память, если вы подумаете о том, чтобы сделать передачу напрямую из сгенерированных данных. Кроме этого, у него есть только четыре дополнительные переменные типа float (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

Отказ от ответственности: я разработал это на бумаге, но еще не имел возможности проверить это "in silico" ...

Вы можете избежать циклов или использования массивов ячеек, выполнив некоторые хитрые кумулятивные суммы, проиндексировав и рассчитав средние и стандартные отклонения самостоятельно. Вот некоторый код, который, я считаю, будет работать, хотя я не уверен, как он сочетается по скорости с другими решениями:

[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 5

Вы можете вычислить indices одновременно, используя bsxfun:

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

Это быстрее, чем зацикливание, но требует их сохранения сразу (компромисс времени и пространства) ..

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);