MATLAB:calcular a média de cada intervalo de 1 minuto de uma série temporal

Question

Eu tenho um monte de séries temporais, cada uma descrita por dois componentes, um vetor de carimbo de data/hora (em segundos) e um vetor de valores medidos.O vetor tempo não é uniforme (ou seja,amostrados em intervalos não regulares)

Estou tentando calcular a média/DP de cada intervalo de valores de 1 minuto (pegue o intervalo de X minutos, calcule sua média, pegue o próximo intervalo, ...).

Minha implementação atual usa loops.Esta é uma amostra do que tenho até agora:

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

Gostaria de saber se existe uma solução vetorizada mais rápida.Isso é importante porque tenho um grande número de séries temporais para processar, cada uma por muito mais tempo do que a amostra mostrada acima.

Qualquer ajuda é bem-vinda.

---

Obrigado a todos pelo feedback.

Eu corrigi o jeito t é gerado para estar sempre aumentando monotonicamente (classificado), isso não era realmente um problema.

Além disso, posso não ter afirmado isso claramente, mas minha intenção era ter uma solução para qualquer intervalo em minutos (1 minuto foi apenas um exemplo)

Answer 1

A única solução lógica parece ser ...

OK. Acho engraçado que para mim haja apenas uma solução lógica, mas muitos outros encontram outras soluções. Independentemente disso, a solução parece simples. Dados os vetores x e t, e um conjunto de pontos de interrupção igualmente espaçados 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)) )';

(Observe que eu classifiquei t acima.)

Eu faria isso em três linhas de código totalmente vetorizadas. Primeiro, se os intervalos fossem arbitrários e potencialmente desiguais no espaçamento, eu usaria o HISTC para determinar quais intervalos a série de dados se enquadra. Dado que são uniformes, apenas faça isso:

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

Novamente, se os elementos de t não fossem classificados por serem classificados, eu teria usado min (t) em vez de t (1). Tendo feito isso, use o Accumarray para reduzir os resultados em um desvio médio e padrão.

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

Answer 2

Você pode tentar criar uma matriz de células e aplicar média e std via Cellfun. É ~ 10% mais lento que a sua solução para 900 entradas, mas ~ 10x mais rápido para 90000 entradas.

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

Nota: Minha solução não fornece exatamente os mesmos resultados que os seus, já que você pula alguns valores de tempo no final (1:60:90 é [1,61]) e, como o início do intervalo não é exatamente o mesmo .

Answer 3

Aqui está uma maneira que usa pesquisa binária.É 6 a 10x mais rápido para 9.900 elementos e cerca de 64x vezes mais rápido para 99.900 elementos.Foi difícil obter tempos confiáveis ​​usando apenas 900 elementos, então não tenho certeza de qual é mais rápido nesse tamanho.Ele quase não usa memória extra se você considerar fazer tx diretamente dos dados gerados.Fora isso, ele possui apenas quatro variáveis ​​flutuantes extras (prevind, first, mid e 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;

Ele usa todas as variáveis ​​que você tinha originalmente.Espero que atenda às suas necessidades.É mais rápido porque leva O(log N) para encontrar os índices com pesquisa binária, mas O(N) para encontrá-los da maneira que você estava fazendo.

Answer 4

Você pode calcular indices de uma só vez usando o BSXFUN:

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

Isso é mais rápido que o loop, mas requer armazená -los de uma só vez (tempo vs troca espacial).

Answer 5

Isenção de responsabilidade: eu trabalhei isso no papel, mas ainda não tive a oportunidade de verificar "em silico" ...

Você pode evitar loops ou usar matrizes de células fazendo algumas somas cumulativas complicadas, indexando e calculando os meios e os desvios padrão. Aqui está algum código que acredito que funcione, embora eu não tenha certeza de como ele empilha em termos de velocidade para as outras soluções:

[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

O acima calcula o desvio padrão usando A simplificação da fórmula encontrada nesta página da Wikipedia.

Answer 6

A mesma resposta acima, mas com o intervalo paramétrico (window_size). Problema com os comprimentos do vetor também resolvidos.

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