MATLAB: media di calcolo di ogni intervallo di 1 minuto di una serie temporale

Question

Ho una serie di serie temporali ciascuna descritta da due componenti, un vettore di data/ora (in secondi) e un vettore di valori misurati. Il vettore del tempo non è uniforme (cioè campionato a intervalli non regolari)

Sto provando a calcolare la media/SD di ciascun intervallo di valori di 1 minuto (prendi X minuto di intervallo, calcola la sua media, prendi il prossimo intervallo, ...).

La mia attuale implementazione utilizza cicli. Questo è un esempio di ciò che ho fino ad ora:

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 

Mi chiedo se esiste una soluzione vettoriale più veloce. Questo è importante perché ho un gran numero di serie temporali da elaborare ognuna molto più a lungo del campione mostrato sopra ..

Qualsiasi aiuto è benvenuto.

---

Grazie a tutti per il feedback.

ho corretto il modo in cui viene generato t di essere sempre monotona crescente (ordinato), questo non era davvero un problema ..

Inoltre, non può aver detto questo, ma evidentemente la mia intenzione era quella di avere una soluzione per qualsiasi intervallo in minuti (1 min è stato solo un esempio)

Answer 1

L'unica soluzione logica sembra essere ...

Ok. Trovo divertente che per me ci sia solo una soluzione logica, ma molti altri trovano altre soluzioni. Indipendentemente da ciò, la soluzione sembra semplice. Dati i vettori x e t, e una serie di equidistanziati punti di rottura 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)))'; 

(noti che risolto t sopra.)

voglio farlo in tre linee completamente vettorializzate di codice. . In primo luogo, in caso di rottura erano arbitrari e potenzialmente diseguale spaziatura, userei histc per determinare quali intervalli la serie di dati cade Dato che sono uniformi, solo fare questo:

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

questo caso, se gli elementi di t non erano noti per essere ordinati, avrei usato min (t) invece di t (1). Fatto ciò, utilizzare accumarray per ridurre i risultati in una deviazione media e standard.

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

Answer 2

Si può provare a creare un array di celle e applicare media e std tramite cellfun. È ~ 10% più lento della soluzione per 900 voci, ma ~ 10 volte più veloce per 90000 voci.

[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: la mia soluzione non dà gli stessi risultati esatti come la tua, da quando si salta un paio di valori di tempo alla fine (1:60:90 è [1,61]), e dal momento che l'inizio di l'intervallo non è esattamente lo stesso.

Answer 3

È possibile calcolare indices tutto in una volta usando bsxfun:

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

Questo è più veloce di looping, ma li richiede la memorizzazione in una volta (il tempo vs spazio compromesso) ..

Answer 4

Ecco un modo che usi binary search. È 6-10 volte più veloce per 9900 elementi e circa 64 volte più veloce per 99900 elementi. È stato difficile ottenere tempi affidabili utilizzando solo 900 elementi, quindi non sono sicuro di quale sia più veloce a quella dimensione. Non utilizza quasi memoria extra se si considera di fare tx direttamente dai dati generati. Oltre a questo ha solo quattro variabili float extra (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; 

Utilizza tutte le variabili che avevate in origine. Spero che sia adatto alle tue esigenze. È più veloce perché richiede O (log N) per trovare gli indici con la ricerca binaria, ma O (N) per trovarli nel modo in cui lo stavi facendo.

Answer 5

Disclaimer: questo ha lavorato su carta, ma non hanno ancora avuto la possibilità di controllare "in silico" ...

Si può essere in grado di evitare loop o utilizzando gli array di celle facendo alcune ingombranti somme cumulative, l'indicizzazione e il calcolo delle medie e delle deviazioni standard.Ecco un po 'di codice che credo funzionerà, anche se non sono sicuro come impila in su velocità-saggio per le altre soluzioni:

[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 

È possibile che questo calcola la deviazione standard utilizzando the simplification of the formula found on this Wikipedia page.

Answer 6

La stessa risposta come sopra ma con l'intervallo parametrico (window_size). Problema risolto anche con le lunghezze dei vettori.

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