Pytanie MATLAB: oblicza średnią każdego 1-minutowego przedziału czasowego


Mam szereg serii czasowych, każdy opisany przez dwa składniki, wektor znacznika czasu (w sekundach) i wektor zmierzonych wartości. Wektor czasu jest niejednolity (tj. Próbkowany w nieregularnych odstępach czasu)

Próbuję obliczyć średnią / SD każdego 1-minutowego przedziału wartości (weź X minut interwału, oblicz jego średnią, weź następny interwał, ...).

Moja obecna implementacja używa pętli. Oto próbka tego, co do tej pory:

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

Zastanawiam się, czy istnieje szybsze wektoryzowane rozwiązanie. Jest to ważne, ponieważ mam dużą liczbę szeregów czasowych do przetworzenia każdego znacznie dłużej niż przykład pokazany powyżej.

Każda pomoc jest mile widziana.


Dziękuję wszystkim za opinie.

Poprawiłem drogę t jest generowany, aby zawsze monotonicznie zwiększać (sortować), to nie był problem.

Być może nie powiedziałem tego wyraźnie, ale moim zamiarem było znalezienie rozwiązania dla dowolnej długości interwału w minutach (1-minuta była tylko przykładem)


10
2018-02-24 01:19


pochodzenie




Odpowiedzi:


Jedynym logicznym rozwiązaniem wydaje się być ...

Ok. Uważam za zabawne, że dla mnie istnieje tylko jedno logiczne rozwiązanie, ale wiele innych znajduje inne rozwiązania. Niezależnie od tego rozwiązanie wydaje się proste. Biorąc pod uwagę wektory x i t oraz zbiór równo rozłożonych punktów przerwania 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)) )';

(Zauważ, że posortowałem t powyżej.)

Zrobiłbym to w trzech całkowicie wektoryzowanych liniach kodu. Po pierwsze, jeśli przerwy byłyby arbitralne i potencjalnie nierówne w rozstawieniu, użyłbym histc do określenia, które interwały ma seria danych. Biorąc pod uwagę, że są one jednorodne, po prostu wykonaj to:

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

Ponownie, gdyby elementy t nie były znane, użyłbym min (t) zamiast t (1). Po wykonaniu tej czynności użyj accumarray, aby zmniejszyć wyniki do średniej i odchylenia standardowego.

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

10
2018-02-24 11:17



+1: Z jakiegoś powodu całkowicie przeoczyłem ACCUMARRAY. - gnovice
dzięki, to jest zarówno zwięzłe i łatwe do odczytania - merv
Nawet nie wiedziałem o akumarobie. Dzięki za wykazanie, jak przydatne może być! - Jonas


Możesz spróbować utworzyć tablicę komórek i zastosować średnią i std za pomocą cellfun. Jest o ~ 10% wolniejsze niż twoje rozwiązanie dla 900 wpisów, ale ~ 10x szybciej dla 90000 wpisów.

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

Uwaga: moje rozwiązanie nie daje dokładnie takich samych wyników jak twoje, ponieważ pomijasz kilka wartości czasu na końcu (1:60:90 to [1,61]), a ponieważ początek przedziału nie jest dokładnie taki sam .


4
2018-02-24 02:25



Dzięki! Mam kilka punktów: [1] masz rację co do sposobu, w jaki generowałem t nie zawsze może być monotonicznie rosnący, nie było to zamierzone! [2] Mimo że wciąż rozszyfrowuję kod, naprawdę potrzebuję, aby długość interwału była sparametryzowana (5 minut to to, nad czym teraz pracuję, ale to powinno być łatwe do zmienienia) ... - merv
[3] prawda jest po obliczeniu stepIdx Trochę się zgubiłem :) Mogłem wyjaśnić co nIdx reprezentuje? Dostaję tę część, w której obliczasz minutową część każdego znacznika czasu, a następnie biorę różnice, aby znaleźć miejsce, w którym się zmienia, wskazując następny 1-minutowy interwał, ale nie mogłem po tym pójść ... - merv
nDX to liczba wyświetleń każdego indeksu. Potrzebuję tego, aby móc używać mat2cell, który dystrybuuje pierwsze n wartości do pierwszej komórki, drugie n wartości w drugiej komórce itd., Grupując w ten sposób indeksy, które należą do każdego przedziału czasu. Mam nadzieję, że dodatkowe komentarze pomogą w zrozumieniu. Przepraszamy za pisanie trudnego do odczytania kodu. Powinienem (pracowałem) nad czymś innym, więc odpowiedziałem na to w pośpiechu :) - Jonas
dziękuję, doceniam pomoc. Co, jeśli chcę mieć inną długość interwału (nie tylko 1 minutę)? kolejne podstępne sztuczki nie będą działać, jakikolwiek pomysł, jak zmienić to na dowolną X-minutę? - merv
Zobacz komentarze w funkcji. Po prostu twórz tidx, dzieląc przez s sekund, gdzie s jest liczbą sekund w przedziale. - Jonas


Oto sposób, który wykorzystuje wyszukiwanie binarne. Jest 6-10-krotnie szybszy dla 9900 elementów i około 64 razy szybszy dla 99900 elementów. Trudno było uzyskać wiarygodne czasy używając tylko 900 elementów, więc nie jestem pewien, który jest szybszy w tym rozmiarze. Nie wykorzystuje prawie żadnej dodatkowej pamięci, jeśli rozważasz wykonanie TX bezpośrednio z wygenerowanych danych. Poza tym ma tylko cztery dodatkowe zmienne zmiennopozycyjne (prevind, first, mid i 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;

Wykorzystuje wszystkie zmienne, które pierwotnie miałeś. Mam nadzieję, że pasuje do twoich potrzeb. Jest szybszy, ponieważ zajmuje O (log N), aby znaleźć indeksy z wyszukiwaniem binarnym, ale O (N), aby znaleźć je w taki sposób, w jaki to robiłeś.


3
2018-02-24 05:40



Powinno to być jeszcze szybsze, jeśli wstępnie zmienisz mu i sd, zamiast rozwijać je wewnątrz pętli. - Jonas
@ Jon: Myślałem, że to będzie implikowane, ponieważ było w kodzie pytającego. Ma to jedynie zastąpić ostatnie 5 linii kodu pytającego. Myślałem, że ostatnie 5 linii to te wolne. - Justin Peel
Czy wyszukiwanie binarne (z pętlami) jest szybsze niż wektor wektorowy, z którym zacząłem? - merv
@merv Tak, porównania czasów, które zrobiłem, zostały porównane z wersją zamieszczoną w pytaniu. Wektoryzacja doprowadza cię tylko do tej pory; przyspiesza to wykonywanie tej samej operacji na każdym elemencie niż w pętli, ale nadal wykonuje operację na każdym elemencie. Oznacza to, że twoja metoda wykonuje 14 * 900 = 12600 porównań (dla 900 elementów), aby znaleźć indeksy, które przyjmą średnią lub std, podczas gdy ta metoda wyszukiwania binarnego wykonuje mniej niż 14 * porównań log_2 (900) = 140 porównań, aby je znaleźć. Różnica się tylko pogarsza wraz ze wzrostem liczebności i rosnącą liczbą interwałów. - Justin Peel


Możesz obliczyć indices wszystkie za pomocą bsxfun:

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

Jest to szybsze niż zapętlenie, ale wymaga przechowywania wszystkich naraz (kompromis czasu i przestrzeni).


2
2018-02-24 04:11



Podoba mi się ten. Jedynym problemem jest to, że nie mogę używać indeksów bezpośrednio bez pętli for: robienie x(indices) nie działa, zamiast tego muszę: for i=1:N, x(indices(:,i)), end - merv


Zastrzeżenie: opracowałem to na papierze, ale nie miałem jeszcze okazji sprawdzić go "in silico" ...

Możesz być w stanie uniknąć pętli lub korzystania z macierzy komórek, wykonując skomplikowane sumy, indeksowanie i obliczając środki i standardowe odchylenia samodzielnie. Oto kod, który, jak sądzę, zadziała, chociaż nie jestem pewny, jak układa się on z szybkością do innych rozwiązań:

[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

Powyższe oblicza standardowe odchylenie za pomocą uproszczenie formuły znajdującej się na tej stronie Wikipedii.


2
2018-02-24 06:31



Dzięki za odpowiedź, myślę, że byłoby interesujące porównanie czasu z innymi rozwiązaniami. - merv


Ta sama odpowiedź co powyżej, ale z interwałem parametrycznym (window_size). Problem z rozwiązaniami wektorowymi.

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

0
2017-12-02 14:37