Matlab: compute moyenne de chaque intervalle de 1 minute d'une série chronologique

Question

I ai un tas de fois chaque série décrite par deux éléments, un vecteur d'horodatage (en secondes), et un vecteur de valeurs mesurées. Le vecteur de temps est non-uniforme (à savoir échantillonné à des intervalles non réguliers)

Je suis en train de calculer la moyenne / SD de chaque intervalle de valeurs 1 minutes (intervalle de prendre X minutes, calculer sa moyenne, prendre l'intervalle, ...).

Mon implémentation 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 pour traiter chaque bien plus que l'échantillon ci-dessus ..

Toute aide est la bienvenue.

---

Merci à tous pour les commentaires.

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

En outre, je ne l'ai dit clairement, mais mon intention était d'avoir une solution pour toute longueur d'intervalle en minutes (1 min était juste un exemple)

Answer 1

La seule solution logique semble être ...

Ok. Je trouve drôle que pour moi il n'y a qu'une seule solution logique, mais beaucoup d'autres à trouver d'autres solutions. Quoiqu'il en soit, la solution ne semble simple. Etant donné les vecteurs x et t, et un ensemble de points de rupture égale distance 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 ferais cela dans 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 ne sont pas connus pour être triés, je l'aurais utilisé min (t) au lieu de t (1). Après avoir fait cela, utilisez accumArray pour réduire les résultats en moyenne et écart-type.

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

Answer 2

Vous pouvez essayer de créer un réseau de cellules et d'appliquer la moyenne et std via cellfun. Il est environ 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 est pas exactement la même chose.

Answer 3

Voici une façon qui utilise recherche binaire . 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 qui est plus rapide à cette taille. Il utilise presque pas de mémoire supplémentaire si vous envisager de faire tx directement à partir des données générées. Autre que celui qu'il a juste quatre variables flottantes supplémentaires (prevind, d'abord, milieu et dernier).

% 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 convient à vos besoins. Il est plus rapide, car il faut O (log N) pour trouver les indices avec la recherche binaire, mais O (N) pour les trouver la façon dont vous le faites.

Answer 4

Vous pouvez calculer indices à la fois en utilisant bsxfun:

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

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

Answer 5

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

Vous pouvez être en mesure d'éviter des boucles ou en utilisant des réseaux cellulaires en faisant des sommes cumulées délicate, l'indexation et le calcul des moyennes et écarts types vous-même. Voici un code que je crois fonctionnera, même si je ne suis pas sûr comment il empile-sage 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 en utilisant la simplification de la formule trouvée sur cette page Wikipedia .

Answer 6

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

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