MATLAB: moyenne de calcul de chaque intervalle de 1 minute d'une série temporelle

Question

J'ai plusieurs séries de temps décrites chacune par deux composantes, un vecteur d'horodatage (en secondes) et un vecteur de valeurs mesurées. Le vecteur temporel est non uniforme (c'est-à-dire échantillonné à intervalles non réguliers)

J'essaie de calculer la moyenne/SD de chaque intervalle de 1 minute de valeurs (prendre X intervalle de minutes, calculer sa moyenne, prendre la prochaine intervalle, ...).

Ma mise en œuvre actuelle utilise des boucles. Ceci est un exemple de ce que j'ai jusqu'à présent:

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 

Je me demande s'il y a une solution plus rapide vectorisé. Ceci est important parce que j'ai un grand nombre de séries chronologiques à traiter chacune beaucoup plus longtemps que l'échantillon montré ci-dessus ..

Toute aide est la bienvenue.

---

Merci à tous pour vos commentaires.

Je corrige la façon dont t est généré pour être toujours de plus en plus de façon monotone (trié), ce n'était pas vraiment un problème ..

En outre, je ne peux pas avoir dit clairement, mais mon intention était d'avoir une solution pour n'importe quelle longueur d'intervalle en minutes (1 minute était juste un exemple)

Answer 1

La seule solution logique semble être ...

Ok. Je trouve amusant que pour moi il n'y ait qu'une seule solution logique, mais que beaucoup d'autres trouvent d'autres solutions. Peu importe, la solution semble simple. Compte tenu des vecteurs x et t, et un ensemble de points de rupture équidistants 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)))'; 

(Notez que je triai t ci-dessus.)

Je ferai en trois lignes entièrement vectorisées de code. . Tout d'abord, si les pauses étaient en espacement arbitraire et potentiellement inégale, je voudrais utiliser histc pour déterminer quels Intervalles les séries de données tombe Étant donné qu'elles sont uniformes, faire ceci:

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

Encore une fois, si les éléments de t n'étaient pas connus pour être triés, j'aurais utilisé min (t) au lieu de t (1). Ayant fait cela, utilisez accumarray pour réduire les résultats dans un écart moyen et standard.

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

Answer 2

Vous pourriez essayer de créer un tableau de cellules et d'appliquer des moyennes et des std via cellfun. C'est ~ 10% plus lent que votre solution pour 900 entrées, mais ~ 10x plus rapide pour 90000 entrées.

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

Note: ma solution ne donne pas les mêmes résultats que la vôtre, puisque vous sauter quelques valeurs de temps à la fin (1:60:90 est [1,61]), et depuis le début de l'intervalle n'est pas exactement le même.

Answer 3

Vous pouvez calculer indices à la fois en utilisant bsxfun:

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

C'est plus rapide que boucle mais nécessite de les stocker à la fois (temps vs compromis entre l'espace) ..

Answer 4

Voici une façon que les utilisations binary search. Il est 6-10x plus rapide pour 9900 éléments et environ 64x fois plus rapide pour 99900 éléments. Il était difficile d'obtenir des temps fiables en utilisant seulement 900 éléments donc je ne suis pas sûr de ce qui est le plus rapide à cette taille. Il n'utilise presque pas de mémoire supplémentaire si vous envisagez de créer tx directement à partir des données générées. Autre que cela, il a juste quatre variables float supplémentaires (prevind, first, mid et 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; 

Il utilise toutes les variables que vous aviez à l'origine. J'espère que cela correspond à vos besoins. C'est plus rapide car il faut O (log N) pour trouver les index avec la recherche binaire, mais O (N) pour les trouver comme vous le faisiez.

Answer 5

Disclaimer: Je travaille sur ce papier, mais n'a pas encore eu l'occasion de le vérifier « in silico » ...

Vous pourriez être en mesure d'éviter les boucles ou en utilisant des réseaux cellulaires en faisant des sommes cumulatives, indexées et calculatrices des moyennes et des écarts-types.Voici un code que je crois fonctionnera, même si je ne suis pas sûr comment il empile terme de vitesse aux autres solutions:

[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 

Le calcule l'écart-type ci-dessus à l'aide the simplification of the formula found on this Wikipedia page.

Answer 6

La même réponse que ci-dessus mais avec l'intervalle paramétrique (window_size). Problème avec les longueurs de vecteurs résolues également.

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