MATLAB: calcular la media de cada intervalo de 1 minuto de una serie de tiempo

Question

Tengo un montón de series de tiempo, cada una descrita por dos componentes, un vector de marca de tiempo (en segundos) y un vector de valores medidos. El vector de tiempo no es uniforme (es decir, muestreado en intervalos no regulares)

Estoy tratando de calcular la media/desviación estándar de cada intervalo de 1 minuto de valores (tome el intervalo de X minutos, calcule su media, tome la siguiente intervalo, ...).

Mi implementación actual utiliza bucles. Esta es una muestra de lo que tengo hasta ahora:

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 

Me pregunto si existe una solución vectorizada más rápida. Esto es importante porque tengo un gran número de series de tiempo para procesar cada una mucho más tiempo que la muestra que se muestra arriba.

Cualquier ayuda es bienvenida.

---

Gracias a todos por los comentarios.

corregí la forma t se genera siempre ser monótona creciente (ordenados), esto no era realmente un problema ..

Además, puede que no haya declarado esto con claridad, pero mi intención era tener una solución para cualquier longitud de intervalo en minutos (1 min fue solo un ejemplo)

Answer 1

La única solución lógica parece ser ...

Ok. Me parece gracioso que para mí haya una sola solución lógica, pero muchos otros encuentran otras soluciones. De todos modos, la solución parece simple. Dados los vectores x y t, y un conjunto de igualmente espaciados tt puntos de ruptura,

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

(Tenga en cuenta que solucionaron t anteriormente.)

Me gustaría hacer esto en tres líneas totalmente vectorizado de código. . En primer lugar, si se rompe la fueron arbitrarias y potencialmente desigual en el espaciamiento, usaría histc para determinar qué intervalos de la serie de datos cae en Dado que son uniformes, simplemente hacer esto:

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

De nuevo, si los elementos de t no se sabía que se clasificaran, habría usado min (t) en lugar de t (1). Una vez hecho esto, use accumarray para reducir los resultados en una media y una desviación estándar.

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

Answer 2

Puede intentar crear una matriz de celdas y aplicar mean y std a través de cellfun. Es ~ 10% más lento que su solución para 900 entradas, pero ~ 10 veces más 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: mi solución no da los mismos resultados que la suya, ya que se salta unos valores de tiempo al final (1:60:90 es [1,61]), y desde el inicio de el intervalo no es exactamente el mismo.

Answer 3

Puede calcular indices todos a la vez usando bsxfun:

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

Esto es más rápido que un bucle, pero todos ellos requiere almacenar a la vez (tiempo vs compensación espacio) ..

Answer 4

Aquí está una manera que los usos binary search. Es 6-10 veces más rápido para 9900 elementos y aproximadamente 64 veces más rápido para 99900 elementos. Fue difícil obtener tiempos confiables usando solo 900 elementos, así que no estoy seguro de cuál es más rápido en ese tamaño. Casi no utiliza memoria extra si considera realizar tx directamente a partir de los datos generados. Aparte de eso, solo tiene cuatro variables flotantes adicionales (anterior, primero, medio y último).

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

Utiliza todas las variables que tenía originalmente. Espero que se adapte a tus necesidades. Es más rápido porque toma O (log N) para encontrar los índices con búsqueda binaria, pero O (N) para encontrarlos de la manera que lo estaba haciendo.

Answer 5

Disclaimer: este trabajado en papel, pero que aún no han tenido la oportunidad de comprobarlo "in silico" ...

Usted puede ser capaz de evitar bucles o la utilización de matrices de células haciendo algunas sumas acumulativas complicadas, indexación y calcular los medios y las desviaciones estándar usted mismo.Aquí hay un código que creo que va a funcionar, aunque estoy seguro de cómo esta se acumula en la velocidad en cuanto a las otras soluciones:

[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 

Lo anterior calcula la desviación estándar utilizando the simplification of the formula found on this Wikipedia page.

Answer 6

la misma respuesta que el anterior pero con el intervalo paramétrico (window_size). Problema con las longitudes de vector resueltas también.

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