Wykrywanie sygnału szczytowego w danych szeregów czasowych w czasie rzeczywistym

242

Aktualizacja: Jak dotąd najlepiej działającym algorytmem jest ten .


To pytanie bada niezawodne algorytmy do wykrywania nagłych szczytów w danych szeregów czasowych w czasie rzeczywistym.

Rozważ następujący zestaw danych:

p = [1 1 1.1 1 0.9 1 1 1.1 1 0.9 1 1.1 1 1 0.9 1 1 1.1 1 1 1 1 1.1 0.9 1 1.1 1 1 0.9 1, ...
     1.1 1 1 1.1 1 0.8 0.9 1 1.2 0.9 1 1 1.1 1.2 1 1.5 1 3 2 5 3 2 1 1 1 0.9 1 1 3, ... 
     2.6 4 3 3.2 2 1 1 0.8 4 4 2 2.5 1 1 1];

(Format Matlab, ale nie chodzi o język, ale o algorytm)

Wykres danych

Wyraźnie widać, że istnieją trzy duże szczyty i niektóre małe szczyty. Ten zestaw danych jest konkretnym przykładem klasy zestawów danych timeseries, o które chodzi w pytaniu. Ta klasa zestawów danych ma dwie ogólne cechy:

  1. Występuje podstawowy hałas o ogólnym znaczeniu
  2. Istnieją duże „ szczyty ” lub „ wyższe punkty danych ”, które znacznie odbiegają od hałasu.

Załóżmy również, że:

  • szerokość pików nie może być wcześniej ustalona
  • wysokość pików wyraźnie i znacząco odbiega od innych wartości
  • używany algorytm musi obliczać w czasie rzeczywistym (więc zmieniaj z każdym nowym punktem danych)

W takiej sytuacji należy skonstruować wartość graniczną, która wyzwala sygnały. Jednak wartość granicy nie może być statyczna i musi być ustalana w czasie rzeczywistym na podstawie algorytmu.


Moje pytanie: jaki jest dobry algorytm do obliczania takich progów w czasie rzeczywistym? Czy istnieją specjalne algorytmy dla takich sytuacji? Jakie są najbardziej znane algorytmy?


Wysoko cenione są solidne algorytmy lub przydatne informacje. (może odpowiedzieć w dowolnym języku: chodzi o algorytm)

Jean-Paul
źródło
5
Tam musi być jakiś absolutny wymóg wysokość za to, że szczyt w uzupełnieniu do wymogów już podane. W przeciwnym razie pik w czasie 13 należy uznać za pik. (Równoważnie: jeśli w przyszłości szczyty
wzrosną
Zgadzam się. Załóżmy, że te szczyty to tylko te, które należy wziąć pod uwagę.
Jean-Paul
Być może zadajesz złe pytanie. Zamiast pytać, jak można wykryć bez zwłoki, możesz zapytać, czy możliwe jest niezwłoczne wykrycie określonego rodzaju sygnału, biorąc pod uwagę tylko to, co wiadomo przed tym czasem, lub co należy wiedzieć o sygnale, aby wykryć coś z danym opóźnienie.
hotpaw2
2
Robiłem to, aby wykryć nagłą zmianę natężenia światła na fotosensorze. Zrobiłem to, zmieniając średnią i ignorując wszelkie punkty danych większe niż próg. Należy zauważyć, że ten próg różni się od progu określającego pik. Załóżmy, że uwzględnisz tylko punkty danych, które mieszczą się w granicach jednego stddev do średniej ruchomej, i rozważ te punkty danych z więcej niż trzema stddev jako wartości szczytowe. Ten algorytm sprawdził się bardzo dobrze w naszym kontekście aplikacji.
justhalf
1
O, rozumiem. Nie spodziewałem się tego w formie kodu. Gdybym widział to pytanie wcześniej, prawdopodobnie uzyskałbyś tę odpowiedź znacznie szybciej = D. W każdym razie moją aplikacją było wykrycie, czy fotosensor jest zasłonięty przez źródło światła otoczenia (dlatego potrzebujemy średniej ruchomej, ponieważ źródło światła otoczenia może zmieniać się stopniowo z czasem). Stworzyliśmy to jako grę, w której powinieneś unosić dłoń nad czujnikami według określonego wzoru. = D
justhalf

Odpowiedzi:

334

Solidny algorytm wykrywania pików (z wykorzystaniem wyników Z)

Wymyśliłem algorytm, który działa bardzo dobrze dla tego rodzaju zestawów danych. Opiera się on na zasadzie dyspersji : jeśli nowy punkt danych jest daną x liczbą standardowych odchyleń od pewnej średniej ruchomej, algorytm sygnalizuje (zwany również z-score ). Algorytm jest bardzo solidny, ponieważ konstruuje oddzielną średnią ruchomą i odchylenie, tak aby sygnały nie uszkodziły progu. Przyszłe sygnały są zatem identyfikowane z mniej więcej taką samą dokładnością, niezależnie od ilości wcześniejszych sygnałów. Algorytm bierze 3 wejścia: lag = the lag of the moving window, threshold = the z-score at which the algorithm signalsi influence = the influence (between 0 and 1) of new signals on the mean and standard deviation. Na przykład a lagz 5 użyje ostatnich 5 obserwacji do wygładzenia danych. ZAthresholdz 3,5 zasygnalizuje, że punkt danych znajduje się 3,5 standardowego odchylenia od średniej ruchomej. A influencez 0,5 daje sygnałowi połowę wpływu, jaki mają normalne punkty danych. Podobnie, influencez 0 całkowicie ignoruje sygnały do ​​ponownego obliczenia nowego progu. Wpływ 0 jest zatem najbardziej niezawodną opcją (ale zakłada stacjonarność ); ustawienie opcji wpływu na 1 jest najmniej niezawodne. W przypadku danych niestacjonarnych opcję wpływu należy zatem umieścić gdzieś pomiędzy 0 a 1.

Działa w następujący sposób:

Pseudo kod

# Let y be a vector of timeseries data of at least length lag+2
# Let mean() be a function that calculates the mean
# Let std() be a function that calculates the standard deviaton
# Let absolute() be the absolute value function

# Settings (the ones below are examples: choose what is best for your data)
set lag to 5;          # lag 5 for the smoothing functions
set threshold to 3.5;  # 3.5 standard deviations for signal
set influence to 0.5;  # between 0 and 1, where 1 is normal influence, 0.5 is half

# Initialize variables
set signals to vector 0,...,0 of length of y;   # Initialize signal results
set filteredY to y(1),...,y(lag)                # Initialize filtered series
set avgFilter to null;                          # Initialize average filter
set stdFilter to null;                          # Initialize std. filter
set avgFilter(lag) to mean(y(1),...,y(lag));    # Initialize first value
set stdFilter(lag) to std(y(1),...,y(lag));     # Initialize first value

for i=lag+1,...,t do
  if absolute(y(i) - avgFilter(i-1)) > threshold*stdFilter(i-1) then
    if y(i) > avgFilter(i-1) then
      set signals(i) to +1;                     # Positive signal
    else
      set signals(i) to -1;                     # Negative signal
    end
    # Reduce influence of signal
    set filteredY(i) to influence*y(i) + (1-influence)*filteredY(i-1);
  else
    set signals(i) to 0;                        # No signal
    set filteredY(i) to y(i);
  end
  # Adjust the filters
  set avgFilter(i) to mean(filteredY(i-lag),...,filteredY(i));
  set stdFilter(i) to std(filteredY(i-lag),...,filteredY(i));
end

Poniżej znajdują się praktyczne zasady wyboru dobrych parametrów dla danych.


Próbny

Demonstracja solidnego algorytmu progowego

Kod Matlab dla tego dema można znaleźć tutaj . Aby skorzystać z wersji demo, po prostu uruchom ją i utwórz szereg czasowy, klikając górny wykres. Algorytm zaczyna działać po narysowaniu lagliczby obserwacji.


Wynik

W przypadku pierwotnego pytania ten algorytm daje następujące dane wyjściowe przy użyciu następujących ustawień lag = 30, threshold = 5, influence = 0:

Przykład algorytmu progowego


Implementacje w różnych językach programowania:


Praktyczne zasady konfigurowania algorytmu

lag: parametr lag określa, o ile dane zostaną wygładzone i jak adaptacyjny jest algorytm do zmian średniej długoterminowej danych. Im bardziej stacjonarne są twoje dane, tym więcej opóźnień powinieneś uwzględnić (powinno to poprawić niezawodność algorytmu). Jeśli Twoje dane zawierają trendy zmieniające się w czasie, powinieneś rozważyć, jak szybko algorytm musi się dostosować do tych trendów. To znaczy, jeśli ustawisz lagna 10, zajmie to 10 „okresów”, zanim próg algorytmu zostanie dostosowany do wszelkich systematycznych zmian średniej długoterminowej. Wybierz więc lagparametr na podstawie trendu zachowania danych i tego, jak elastyczny ma być algorytm.

influence: ten parametr określa wpływ sygnałów na próg detekcji algorytmu. Jeśli ustawione na 0, sygnały nie mają wpływu na wartość progową, tak że przyszłe sygnały są wykrywane na podstawie wartości progowej, która jest obliczana ze średniej i odchylenia standardowego, na które nie mają wpływu sygnały z przeszłości. Innym sposobem myślenia o tym jest to, że jeśli ustawisz wpływ na 0, domyślnie przyjmiesz stacjonarność (tj. Bez względu na liczbę sygnałów, szeregi czasowe zawsze wracają do tej samej średniej w dłuższej perspektywie). Jeśli tak nie jest, parametr wpływu należy umieścić gdzieś pomiędzy 0 a 1, w zależności od stopnia, w jakim sygnały mogą systematycznie wpływać na zmieniający się w czasie trend danych. Np. Jeśli sygnały prowadzą do pęknięcia strukturalnego średniej długoterminowej szeregu czasowego parametr wpływu należy ustawić na wysokim poziomie (blisko 1), aby próg mógł szybko dostosować się do tych zmian.

threshold: parametr progowy to liczba odchyleń standardowych od średniej ruchomej, powyżej której algorytm sklasyfikuje nowy punkt danych jako sygnał. Na przykład, jeśli nowy punkt danych ma 4,0 odchylenia standardowe powyżej średniej ruchomej, a parametr progowy ustawiony jest na 3,5, algorytm zidentyfikuje punkt danych jako sygnał. Ten parametr należy ustawić na podstawie oczekiwanej liczby sygnałów. Na przykład, jeśli dane są normalnie dystrybuowane, próg (lub: wynik-Z) wynoszący 3,5 odpowiada prawdopodobieństwu sygnalizacyjnemu wynoszącemu 0,00047 (z tej tabeli), co oznacza, że ​​oczekujesz sygnału raz na 2128 punktów danych (1 / 0,00047). Próg wpływa zatem bezpośrednio na czułość algorytmu, a tym samym także na częstotliwość sygnałów. Zbadaj własne dane i określ rozsądny próg, który sygnalizuje algorytmowi, kiedy tego chcesz (może być potrzebna próba i błąd, aby uzyskać dobry próg dla twojego celu).


OSTRZEŻENIE: Powyższy kod zawsze zapętla wszystkie punkty danych przy każdym uruchomieniu. Wdrażając ten kod, pamiętaj o podzieleniu obliczeń sygnału na osobną funkcję (bez pętli). Wtedy, gdy nowy Datapoint przybywa, aktualizacja filteredY, avgFiltera stdFilterraz. Nie przeliczaj sygnałów dla wszystkich danych za każdym razem, gdy pojawia się nowy punkt danych (jak w powyższym przykładzie), który byłby wyjątkowo nieefektywny i powolny!

Inne sposoby modyfikacji algorytmu (w celu potencjalnych ulepszeń) to:

  1. Użyj średniej zamiast średniej
  2. Zastosuj solidną miarę skali , taką jak MAD, zamiast standardowego odchylenia
  3. Użyj marginesu sygnalizacyjnego, aby sygnał nie zmieniał się zbyt często
  4. Zmień sposób działania parametru wpływu
  5. Różnie traktuj sygnały w górę i w dół (leczenie asymetryczne)
  6. Utwórz osobny influenceparametr dla średniej i standardowej wartości ( tak jak w tym tłumaczeniu Swift )

(Znane) cytaty naukowe do tej odpowiedzi StackOverflow:

Inna praca przy użyciu algorytmu

Inne zastosowania tego algorytmu

Linki do innych algorytmów wykrywania pików


Jeśli gdzieś skorzystasz z tej funkcji, proszę o kredyt lub tę odpowiedź. Jeśli masz jakieś pytania dotyczące tego algorytmu, opublikuj je w komentarzach poniżej lub skontaktuj się ze mną na LinkedIn .


Jean-Paul
źródło
Link do movingstd jest zepsuty, ale jego opis można znaleźć tutaj
Phylliida
@reasra Okazuje się, że funkcja nie potrzebuje ruchomego odchylenia standardowego po przepisaniu. Można go teraz używać z prostymi, wbudowanymi funkcjami Matlaba :)
Jean-Paul,
1
Próbuję kodu Matlab dla niektórych danych akcelerometru, ale z jakiegoś powodu thresholdwykres staje się płaską zieloną linią po dużym skoku do 20 w danych i pozostaje taki przez resztę wykresu ... Jeśli Usuwam sike, tak się nie dzieje, więc wydaje się, że jest to spowodowane skokiem danych. Masz pomysł, co się dzieje? Jestem nowicjuszem w Matlabie, więc nie mogę tego rozgryźć ...
Magnus W
@BadCash Czy możesz podać przykład (z danymi)? Być może zadaj swoje własne pytanie tutaj na SO i powiedz mi link?
Jean-Paul,
2
Istnieje wiele sposobów na ulepszenie tego algo, więc bądź kreatywny (inne leczenie w górę / w dół; mediana zamiast średniej; solidny standard; pisanie kodu jako funkcji wydajnej pamięci; margines progowy, aby sygnał nie przełączał się zbyt często itp. .).
Jean-Paul,
41

Oto implementacja Python/ numpywygładzonego algorytmu Z-score (patrz odpowiedź powyżej ). Można znaleźć sedno tutaj .

#!/usr/bin/env python
# Implementation of algorithm from https://stackoverflow.com/a/22640362/6029703
import numpy as np
import pylab

def thresholding_algo(y, lag, threshold, influence):
    signals = np.zeros(len(y))
    filteredY = np.array(y)
    avgFilter = [0]*len(y)
    stdFilter = [0]*len(y)
    avgFilter[lag - 1] = np.mean(y[0:lag])
    stdFilter[lag - 1] = np.std(y[0:lag])
    for i in range(lag, len(y)):
        if abs(y[i] - avgFilter[i-1]) > threshold * stdFilter [i-1]:
            if y[i] > avgFilter[i-1]:
                signals[i] = 1
            else:
                signals[i] = -1

            filteredY[i] = influence * y[i] + (1 - influence) * filteredY[i-1]
            avgFilter[i] = np.mean(filteredY[(i-lag+1):i+1])
            stdFilter[i] = np.std(filteredY[(i-lag+1):i+1])
        else:
            signals[i] = 0
            filteredY[i] = y[i]
            avgFilter[i] = np.mean(filteredY[(i-lag+1):i+1])
            stdFilter[i] = np.std(filteredY[(i-lag+1):i+1])

    return dict(signals = np.asarray(signals),
                avgFilter = np.asarray(avgFilter),
                stdFilter = np.asarray(stdFilter))

Poniżej znajduje się test na tym samym zestawie danych, który daje taki sam wykres jak w oryginalnej odpowiedzi na R/Matlab

# Data
y = np.array([1,1,1.1,1,0.9,1,1,1.1,1,0.9,1,1.1,1,1,0.9,1,1,1.1,1,1,1,1,1.1,0.9,1,1.1,1,1,0.9,
       1,1.1,1,1,1.1,1,0.8,0.9,1,1.2,0.9,1,1,1.1,1.2,1,1.5,1,3,2,5,3,2,1,1,1,0.9,1,1,3,
       2.6,4,3,3.2,2,1,1,0.8,4,4,2,2.5,1,1,1])

# Settings: lag = 30, threshold = 5, influence = 0
lag = 30
threshold = 5
influence = 0

# Run algo with settings from above
result = thresholding_algo(y, lag=lag, threshold=threshold, influence=influence)

# Plot result
pylab.subplot(211)
pylab.plot(np.arange(1, len(y)+1), y)

pylab.plot(np.arange(1, len(y)+1),
           result["avgFilter"], color="cyan", lw=2)

pylab.plot(np.arange(1, len(y)+1),
           result["avgFilter"] + threshold * result["stdFilter"], color="green", lw=2)

pylab.plot(np.arange(1, len(y)+1),
           result["avgFilter"] - threshold * result["stdFilter"], color="green", lw=2)

pylab.subplot(212)
pylab.step(np.arange(1, len(y)+1), result["signals"], color="red", lw=2)
pylab.ylim(-1.5, 1.5)
pylab.show()
R Kiselev
źródło
Tutaj „y” to tak naprawdę sygnał, a „sygnały” to zbiór punktów danych, czy mam rację, rozumiejąc?
TheTank
1
@TheTank yoznacza macierz danych można przejść w, signalsto +1albo -1tablica wyjście, które wskazują na każdym Datapoint y[i]czy Datapoint jest „znacząca szczytowa” Biorąc pod uwagę ustawienia użyć.
Jean-Paul,
23

Jednym z podejść jest wykrywanie pików na podstawie następującej obserwacji:

  • Czas t jest szczytem, ​​jeśli (y (t)> y (t-1)) i& (y (t)> y (t + 1))

Unika fałszywych trafień, czekając, aż trend wzrostowy się skończy. Nie jest to dokładnie „w czasie rzeczywistym” w tym sensie, że spóźni się o szczyt o jeden dt. czułość można kontrolować, wymagając marginesu dla porównania. Istnieje kompromis między głośnym wykrywaniem a opóźnieniem czasowym wykrycia. Możesz wzbogacić model, dodając więcej parametrów:

  • pik jeśli (y (t) - y (t-dt)> m) i& (y (t) - y (t + dt)> m)

gdzie DT i m są parametrami do kontroli czułości w funkcji czasu opóźnienia

Oto, co otrzymujesz dzięki wspomnianemu algorytmowi: wprowadź opis zdjęcia tutaj

oto kod do odtworzenia wykresu w Pythonie:

import numpy as np
import matplotlib.pyplot as plt
input = np.array([ 1. ,  1. ,  1. ,  1. ,  1. ,  1. ,  1. ,  1.1,  1. ,  0.8,  0.9,
    1. ,  1.2,  0.9,  1. ,  1. ,  1.1,  1.2,  1. ,  1.5,  1. ,  3. ,
    2. ,  5. ,  3. ,  2. ,  1. ,  1. ,  1. ,  0.9,  1. ,  1. ,  3. ,
    2.6,  4. ,  3. ,  3.2,  2. ,  1. ,  1. ,  1. ,  1. ,  1. ])
signal = (input > np.roll(input,1)) & (input > np.roll(input,-1))
plt.plot(input)
plt.plot(signal.nonzero()[0], input[signal], 'ro')
plt.show()

Ustawiając m = 0.5, możesz uzyskać czystszy sygnał z tylko jednym fałszywie dodatnim: wprowadź opis zdjęcia tutaj

Aha
źródło
Wcześniej = lepsze, więc wszystkie piki są znaczące. Dzięki! Bardzo fajny!
Jean-Paul
Jak mogę zmienić czułość?
Jean-Paul
Mogę wymyślić dwa podejścia: 1: ustaw m na większą wartość, aby wykryć tylko większe piki. 2: zamiast obliczać y (t) - y (t-dt) (i y (t) - y (t + dt)), całkujesz od t-dt do t (i t do t + dt).
aha
2
Według jakich kryteriów odrzucasz pozostałe 7 szczytów?
hotpaw2
4
Występuje problem z płaskimi pikami, ponieważ to, co robisz, to w zasadzie wykrywanie zbocza 1-D (jak zwijanie sygnału za pomocą [1 0 -1])
ben
18

W przetwarzaniu sygnału wykrywanie pików często odbywa się za pomocą transformaty falkowej. Zasadniczo wykonujesz dyskretną transformatę falkową na danych szeregów czasowych. Zerowane przecięcia w zwracanych współczynnikach szczegółowych będą odpowiadać pikom w sygnale szeregów czasowych. Otrzymujesz różne amplitudy pików wykrywane przy różnych poziomach współczynnika szczegółowości, co daje rozdzielczość wielopoziomową.

Cklin
źródło
1
Twoja odpowiedź pozwoliła mi przejść do tego artykułu i tej odpowiedzi, która pomoże mi skonstruować dobry algorytm dla mojej implementacji. Dzięki!
Jean-Paul
@cklin Czy możesz wyjaśnić, w jaki sposób obliczasz przejścia przez zero falek, ponieważ nie są one w tej samej skali czasowej, co oryginalne szeregi czasowe. Jakieś referencje na temat tego użycia?
horaceT
11

Podjęliśmy próbę zastosowania wygładzonego algorytmu Z-score w naszym zestawie danych, co powoduje albo nadwrażliwość, albo nadwrażliwość (w zależności od tego, jak parametry są dostrojone), z niewielkim pośrednim poziomem. W sygnalizacji ruchu na naszej stronie zaobserwowaliśmy linię bazową niskiej częstotliwości, która reprezentuje cykl dzienny, a nawet przy najlepszych możliwych parametrach (pokazanych poniżej), nadal występowała, szczególnie czwartego dnia, ponieważ większość punktów danych jest rozpoznawana jako anomalia .

Opierając się na oryginalnym algorytmie Z-score, znaleźliśmy sposób na rozwiązanie tego problemu przez filtrowanie wsteczne. Szczegóły zmodyfikowanego algorytmu i jego zastosowania do komercyjnego przypisywania ruchu telewizyjnego są publikowane na blogu naszego zespołu .

wprowadź opis zdjęcia tutaj

jbm
źródło
Fajnie jest zobaczyć, że algorytm był punktem wyjścia dla bardziej zaawansowanej wersji. Twoje dane mają bardzo szczególny wzór, więc rzeczywiście sensowniej byłoby najpierw usunąć wzór za pomocą innej techniki, a następnie zastosować algo do reszt. Alternatywnie, możesz użyć wyśrodkowanego zamiast opóźnionego okna, aby obliczyć średnią / odchylenie standardowe. Kolejny komentarz: twoje rozwiązanie przesuwa się z prawej na lewą, aby zidentyfikować skoki, ale nie jest to możliwe w aplikacjach w czasie rzeczywistym (dlatego oryginalne algo jest tak uproszczone, ponieważ przyszłe informacje są niedostępne).
Jean-Paul
10

W topologii obliczeniowej idea trwałej homologii prowadzi do wydajnego - szybkiego jak sortowanie liczb - rozwiązania. Nie tylko wykrywa szczyty, ale w naturalny sposób określa ilościowo „znaczenie” pików, co pozwala wybrać piki, które są dla Ciebie znaczące.

Podsumowanie algorytmu. W ustawieniu 1-wymiarowym (szeregi czasowe, sygnał wartości rzeczywistej) algorytm można łatwo opisać na poniższym rysunku:

Najbardziej trwałe szczyty

Pomyśl o wykresie funkcji (lub jego zestawie podpoziomowym) jako krajobrazie i rozważ obniżający się poziom wody, poczynając od poziomu nieskończoności (lub 1,8 na tym zdjęciu). Podczas gdy poziom spada, na lokalnych maksymach wyskakują. W lokalnych minimach wyspy te łączą się ze sobą. Jednym ze szczegółów tego pomysłu jest to, że wyspa, która pojawiła się później, połączona jest ze starszą wyspą. „Trwałość” wyspy to czas jej narodzin minus czas jej śmierci. Długości niebieskich słupków obrazują trwałość, która jest wyżej wspomnianym „znaczeniem” piku.

Wydajność. Nietrudno jest znaleźć implementację działającą w czasie liniowym - w rzeczywistości jest to pojedyncza, prosta pętla - po posortowaniu wartości funkcji. Dlatego wdrożenie to powinno być szybkie w praktyce i łatwe do wdrożenia.

Bibliografia. Podsumowanie całej historii i odniesienia do motywacji z trwałej homologii (pole w obliczeniowej topologii algebraicznej) można znaleźć tutaj: https://www.sthu.org/blog/13-perstopology-peakdetection/index.html

S. Huber
źródło
Ten algorytm jest znacznie szybszy i dokładniejszy niż na przykład scipy.signal.find_peaks. Dla „rzeczywistej” serii czasowej z 1053896 punktami danych wykrył 137516 pików (13%). Kolejność pików (najpierw najbardziej znacząca) pozwala na wydobycie najbardziej znaczących pików. Zapewnia początek, szczyt i koniec każdego piku. Działa dobrze z hałaśliwymi danymi.
Vinh
Przez dane w czasie rzeczywistym rozumiesz tak zwany algorytm online, w którym punkty danych są odbierane za każdym razem. Znaczenie piku może zależeć od wartości w przyszłości. Byłoby miło rozszerzyć algorytm, aby stał się on-line, modyfikując wcześniejsze wyniki bez poświęcania zbytniej złożoności czasu.
S. Huber,
9

Znalazłem inny algorytm GH Palshikara w Prostych algorytmach wykrywania pików w szeregach czasowych .

Algorytm wygląda następująco:

algorithm peak1 // one peak detection algorithms that uses peak function S1 

input T = x1, x2, …, xN, N // input time-series of N points 
input k // window size around the peak 
input h // typically 1 <= h <= 3 
output O // set of peaks detected in T 

begin 
O = empty set // initially empty 

    for (i = 1; i < n; i++) do
        // compute peak function value for each of the N points in T 
        a[i] = S1(k,i,xi,T); 
    end for 

    Compute the mean m' and standard deviation s' of all positive values in array a; 

    for (i = 1; i < n; i++) do // remove local peaks which are “small” in global context 
        if (a[i] > 0 && (a[i] – m') >( h * s')) then O = O + {xi}; 
        end if 
    end for 

    Order peaks in O in terms of increasing index in T 

    // retain only one peak out of any set of peaks within distance k of each other 

    for every adjacent pair of peaks xi and xj in O do 
        if |j – i| <= k then remove the smaller value of {xi, xj} from O 
        end if 
    end for 
end

Zalety

  • Artykuł zawiera 5 różnych algorytmów do wykrywania pików
  • Algorytmy działają na nieprzetworzonych danych szeregów czasowych (nie jest wymagane wygładzanie)

Niedogodności

  • Trudne do ustalenia kih wcześniej
  • Piki nie mogą być płaskie (jak trzeci pik w moich danych testowych)

Przykład:

wprowadź opis zdjęcia tutaj

Jean-Paul
źródło
Właściwie interesujący papier. Jego zdaniem S4 wydaje się lepszą funkcją do wykorzystania. Ale ważniejsze jest wyjaśnienie, kiedy k <i <Nk nie jest prawdziwe. Jak zdefiniować funkcję S1 (S2 ... potem po prawej
daniels_pa
8

Oto implementacja algorytmu wygładzonego z-score (powyżej) w Golang. Zakłada kawałek []int16(16-bitowe próbki PCM). Tutaj znajdziesz sedno .

/*
Settings (the ones below are examples: choose what is best for your data)
set lag to 5;          # lag 5 for the smoothing functions
set threshold to 3.5;  # 3.5 standard deviations for signal
set influence to 0.5;  # between 0 and 1, where 1 is normal influence, 0.5 is half
*/

// ZScore on 16bit WAV samples
func ZScore(samples []int16, lag int, threshold float64, influence float64) (signals []int16) {
    //lag := 20
    //threshold := 3.5
    //influence := 0.5

    signals = make([]int16, len(samples))
    filteredY := make([]int16, len(samples))
    for i, sample := range samples[0:lag] {
        filteredY[i] = sample
    }
    avgFilter := make([]int16, len(samples))
    stdFilter := make([]int16, len(samples))

    avgFilter[lag] = Average(samples[0:lag])
    stdFilter[lag] = Std(samples[0:lag])

    for i := lag + 1; i < len(samples); i++ {

        f := float64(samples[i])

        if float64(Abs(samples[i]-avgFilter[i-1])) > threshold*float64(stdFilter[i-1]) {
            if samples[i] > avgFilter[i-1] {
                signals[i] = 1
            } else {
                signals[i] = -1
            }
            filteredY[i] = int16(influence*f + (1-influence)*float64(filteredY[i-1]))
            avgFilter[i] = Average(filteredY[(i - lag):i])
            stdFilter[i] = Std(filteredY[(i - lag):i])
        } else {
            signals[i] = 0
            filteredY[i] = samples[i]
            avgFilter[i] = Average(filteredY[(i - lag):i])
            stdFilter[i] = Std(filteredY[(i - lag):i])
        }
    }

    return
}

// Average a chunk of values
func Average(chunk []int16) (avg int16) {
    var sum int64
    for _, sample := range chunk {
        if sample < 0 {
            sample *= -1
        }
        sum += int64(sample)
    }
    return int16(sum / int64(len(chunk)))
}
Xeoncross
źródło
@ Jean-Paul Nie jestem całkowicie pewien, czy wszystko jest w porządku, więc mogą występować błędy.
Xeoncross,
1
Czy próbowałeś zreplikować przykładowe dane wyjściowe z Matlab / R? To powinno być dobrym potwierdzeniem jakości.
Jean-Paul,
7

Oto implementacja wygładzonego algorytmu z-score C ++ z tej odpowiedzi

std::vector<int> smoothedZScore(std::vector<float> input)
{   
    //lag 5 for the smoothing functions
    int lag = 5;
    //3.5 standard deviations for signal
    float threshold = 3.5;
    //between 0 and 1, where 1 is normal influence, 0.5 is half
    float influence = .5;

    if (input.size() <= lag + 2)
    {
        std::vector<int> emptyVec;
        return emptyVec;
    }

    //Initialise variables
    std::vector<int> signals(input.size(), 0.0);
    std::vector<float> filteredY(input.size(), 0.0);
    std::vector<float> avgFilter(input.size(), 0.0);
    std::vector<float> stdFilter(input.size(), 0.0);
    std::vector<float> subVecStart(input.begin(), input.begin() + lag);
    avgFilter[lag] = mean(subVecStart);
    stdFilter[lag] = stdDev(subVecStart);

    for (size_t i = lag + 1; i < input.size(); i++)
    {
        if (std::abs(input[i] - avgFilter[i - 1]) > threshold * stdFilter[i - 1])
        {
            if (input[i] > avgFilter[i - 1])
            {
                signals[i] = 1; //# Positive signal
            }
            else
            {
                signals[i] = -1; //# Negative signal
            }
            //Make influence lower
            filteredY[i] = influence* input[i] + (1 - influence) * filteredY[i - 1];
        }
        else
        {
            signals[i] = 0; //# No signal
            filteredY[i] = input[i];
        }
        //Adjust the filters
        std::vector<float> subVec(filteredY.begin() + i - lag, filteredY.begin() + i);
        avgFilter[i] = mean(subVec);
        stdFilter[i] = stdDev(subVec);
    }
    return signals;
}
ćwiek
źródło
2
Zastrzeżenie: ta implementacja nie zapewnia metody obliczania średniej i odchylenia standardowego. W przypadku C ++ 11 łatwą metodę można znaleźć tutaj: stackoverflow.com/a/12405793/3250829
rayryeng
6

Ten problem wygląda podobnie do tego, który napotkałem na kursie systemów hybrydowych / wbudowanych, ale był związany z wykrywaniem błędów, gdy sygnał wejściowy z czujnika jest głośny. Zastosowaliśmy filtr Kalmana, aby oszacować / przewidzieć ukryty stan systemu, a następnie zastosowaliśmy analizę statystyczną, aby określić prawdopodobieństwo wystąpienia usterki . Pracowaliśmy z układami liniowymi, ale istnieją warianty nieliniowe. Pamiętam, że to podejście było zaskakująco adaptacyjne, ale wymagało modelu dynamiki systemu.

Peter G.
źródło
Filtr Kalmana jest interesujący, ale nie mogę znaleźć odpowiedniego algorytmu do moich celów. Bardzo doceniam odpowiedź i przyjrzę się niektórym dokumentom z wykrywaniem szczytów, takim jak ten, aby sprawdzić, czy mogę nauczyć się na podstawie dowolnego z algorytmów. Dzięki!
Jean-Paul
6

Implementacja C ++

#include <iostream>
#include <vector>
#include <algorithm>
#include <unordered_map>
#include <cmath>
#include <iterator>
#include <numeric>

using namespace std;

typedef long double ld;
typedef unsigned int uint;
typedef std::vector<ld>::iterator vec_iter_ld;

/**
 * Overriding the ostream operator for pretty printing vectors.
 */
template<typename T>
std::ostream &operator<<(std::ostream &os, std::vector<T> vec) {
    os << "[";
    if (vec.size() != 0) {
        std::copy(vec.begin(), vec.end() - 1, std::ostream_iterator<T>(os, " "));
        os << vec.back();
    }
    os << "]";
    return os;
}

/**
 * This class calculates mean and standard deviation of a subvector.
 * This is basically stats computation of a subvector of a window size qual to "lag".
 */
class VectorStats {
public:
    /**
     * Constructor for VectorStats class.
     *
     * @param start - This is the iterator position of the start of the window,
     * @param end   - This is the iterator position of the end of the window,
     */
    VectorStats(vec_iter_ld start, vec_iter_ld end) {
        this->start = start;
        this->end = end;
        this->compute();
    }

    /**
     * This method calculates the mean and standard deviation using STL function.
     * This is the Two-Pass implementation of the Mean & Variance calculation.
     */
    void compute() {
        ld sum = std::accumulate(start, end, 0.0);
        uint slice_size = std::distance(start, end);
        ld mean = sum / slice_size;
        std::vector<ld> diff(slice_size);
        std::transform(start, end, diff.begin(), [mean](ld x) { return x - mean; });
        ld sq_sum = std::inner_product(diff.begin(), diff.end(), diff.begin(), 0.0);
        ld std_dev = std::sqrt(sq_sum / slice_size);

        this->m1 = mean;
        this->m2 = std_dev;
    }

    ld mean() {
        return m1;
    }

    ld standard_deviation() {
        return m2;
    }

private:
    vec_iter_ld start;
    vec_iter_ld end;
    ld m1;
    ld m2;
};

/**
 * This is the implementation of the Smoothed Z-Score Algorithm.
 * This is direction translation of https://stackoverflow.com/a/22640362/1461896.
 *
 * @param input - input signal
 * @param lag - the lag of the moving window
 * @param threshold - the z-score at which the algorithm signals
 * @param influence - the influence (between 0 and 1) of new signals on the mean and standard deviation
 * @return a hashmap containing the filtered signal and corresponding mean and standard deviation.
 */
unordered_map<string, vector<ld>> z_score_thresholding(vector<ld> input, int lag, ld threshold, ld influence) {
    unordered_map<string, vector<ld>> output;

    uint n = (uint) input.size();
    vector<ld> signals(input.size());
    vector<ld> filtered_input(input.begin(), input.end());
    vector<ld> filtered_mean(input.size());
    vector<ld> filtered_stddev(input.size());

    VectorStats lag_subvector_stats(input.begin(), input.begin() + lag);
    filtered_mean[lag - 1] = lag_subvector_stats.mean();
    filtered_stddev[lag - 1] = lag_subvector_stats.standard_deviation();

    for (int i = lag; i < n; i++) {
        if (abs(input[i] - filtered_mean[i - 1]) > threshold * filtered_stddev[i - 1]) {
            signals[i] = (input[i] > filtered_mean[i - 1]) ? 1.0 : -1.0;
            filtered_input[i] = influence * input[i] + (1 - influence) * filtered_input[i - 1];
        } else {
            signals[i] = 0.0;
            filtered_input[i] = input[i];
        }
        VectorStats lag_subvector_stats(filtered_input.begin() + (i - lag), filtered_input.begin() + i);
        filtered_mean[i] = lag_subvector_stats.mean();
        filtered_stddev[i] = lag_subvector_stats.standard_deviation();
    }

    output["signals"] = signals;
    output["filtered_mean"] = filtered_mean;
    output["filtered_stddev"] = filtered_stddev;

    return output;
};

int main() {
    vector<ld> input = {1.0, 1.0, 1.1, 1.0, 0.9, 1.0, 1.0, 1.1, 1.0, 0.9, 1.0, 1.1, 1.0, 1.0, 0.9, 1.0, 1.0, 1.1, 1.0,
                        1.0, 1.0, 1.0, 1.1, 0.9, 1.0, 1.1, 1.0, 1.0, 0.9, 1.0, 1.1, 1.0, 1.0, 1.1, 1.0, 0.8, 0.9, 1.0,
                        1.2, 0.9, 1.0, 1.0, 1.1, 1.2, 1.0, 1.5, 1.0, 3.0, 2.0, 5.0, 3.0, 2.0, 1.0, 1.0, 1.0, 0.9, 1.0,
                        1.0, 3.0, 2.6, 4.0, 3.0, 3.2, 2.0, 1.0, 1.0, 0.8, 4.0, 4.0, 2.0, 2.5, 1.0, 1.0, 1.0};

    int lag = 30;
    ld threshold = 5.0;
    ld influence = 0.0;
    unordered_map<string, vector<ld>> output = z_score_thresholding(input, lag, threshold, influence);
    cout << output["signals"] << endl;
}
Animesh Pandey
źródło
6

Kontynuując zaproponowane rozwiązanie @ Jean-Paula, zaimplementowałem jego algorytm w C #

public class ZScoreOutput
{
    public List<double> input;
    public List<int> signals;
    public List<double> avgFilter;
    public List<double> filtered_stddev;
}

public static class ZScore
{
    public static ZScoreOutput StartAlgo(List<double> input, int lag, double threshold, double influence)
    {
        // init variables!
        int[] signals = new int[input.Count];
        double[] filteredY = new List<double>(input).ToArray();
        double[] avgFilter = new double[input.Count];
        double[] stdFilter = new double[input.Count];

        var initialWindow = new List<double>(filteredY).Skip(0).Take(lag).ToList();

        avgFilter[lag - 1] = Mean(initialWindow);
        stdFilter[lag - 1] = StdDev(initialWindow);

        for (int i = lag; i < input.Count; i++)
        {
            if (Math.Abs(input[i] - avgFilter[i - 1]) > threshold * stdFilter[i - 1])
            {
                signals[i] = (input[i] > avgFilter[i - 1]) ? 1 : -1;
                filteredY[i] = influence * input[i] + (1 - influence) * filteredY[i - 1];
            }
            else
            {
                signals[i] = 0;
                filteredY[i] = input[i];
            }

            // Update rolling average and deviation
            var slidingWindow = new List<double>(filteredY).Skip(i - lag).Take(lag+1).ToList();

            var tmpMean = Mean(slidingWindow);
            var tmpStdDev = StdDev(slidingWindow);

            avgFilter[i] = Mean(slidingWindow);
            stdFilter[i] = StdDev(slidingWindow);
        }

        // Copy to convenience class 
        var result = new ZScoreOutput();
        result.input = input;
        result.avgFilter       = new List<double>(avgFilter);
        result.signals         = new List<int>(signals);
        result.filtered_stddev = new List<double>(stdFilter);

        return result;
    }

    private static double Mean(List<double> list)
    {
        // Simple helper function! 
        return list.Average();
    }

    private static double StdDev(List<double> values)
    {
        double ret = 0;
        if (values.Count() > 0)
        {
            double avg = values.Average();
            double sum = values.Sum(d => Math.Pow(d - avg, 2));
            ret = Math.Sqrt((sum) / (values.Count() - 1));
        }
        return ret;
    }
}

Przykładowe użycie:

var input = new List<double> {1.0, 1.0, 1.1, 1.0, 0.9, 1.0, 1.0, 1.1, 1.0, 0.9, 1.0,
    1.1, 1.0, 1.0, 0.9, 1.0, 1.0, 1.1, 1.0, 1.0, 1.0, 1.0, 1.1, 0.9, 1.0, 1.1, 1.0, 1.0, 0.9,
    1.0, 1.1, 1.0, 1.0, 1.1, 1.0, 0.8, 0.9, 1.0, 1.2, 0.9, 1.0, 1.0, 1.1, 1.2, 1.0, 1.5, 1.0,
    3.0, 2.0, 5.0, 3.0, 2.0, 1.0, 1.0, 1.0, 0.9, 1.0, 1.0, 3.0, 2.6, 4.0, 3.0, 3.2, 2.0, 1.0,
    1.0, 0.8, 4.0, 4.0, 2.0, 2.5, 1.0, 1.0, 1.0};

int lag = 30;
double threshold = 5.0;
double influence = 0.0;

var output = ZScore.StartAlgo(input, lag, threshold, influence);
Ocean Airdrop
źródło
1
Cześć @ Jean-Paul. Twoje zdrowie. Tak, przetestowałem dane wyjściowe pod kątem wersji R, aby upewnić się, że są zgodne. Jeszcze raz dziękuję za rozwiązanie tego problemu.
Ocean Airdrop,
Cześć, myślę, że jest błąd w tym kodzie, w metodzie StdDev bierzesz wartości.Count () - 1, czy powinno być -1? Myślę, że chciałbyś mieć liczbę elementów i to właśnie otrzymujesz z wartości.Count ().
Viktor
1
Hmm .. Dobre miejsce. Chociaż pierwotnie przeniosłem algorytm do C #, nigdy go nie użyłem. Prawdopodobnie zastąpiłbym tę całą funkcję wywołaniem biblioteki nuget MathNet. „Install-Package MathNet.Numerics” Posiada wstępnie wbudowane funkcje dla PopulationStandardDeviation () i StandardDeviation (); na przykład. var populacjaStdDev = nowa lista <double> (1,2,3,4) .PopulationStandardDeviation (); var sampleStdDev = nowa lista <double> (1,2,3,4) .StandardDeviation ();
Ocean Airdrop
6

Oto implementacja C wygładzonego wyniku Z @ Jean-Paula dla mikrokontrolera Arduino wykorzystywanego do odczytu odczytów akcelerometru i decydowania, czy kierunek uderzenia pochodzi z lewej, czy z prawej strony. Działa to naprawdę dobrze, ponieważ to urządzenie zwraca odbity sygnał. Oto dane wejściowe do tego algorytmu wykrywania pików z urządzenia - pokazujące uderzenie z prawej strony, a następnie uderzenie z lewej. Możesz zobaczyć początkowy skok, a następnie oscylację czujnika.

wprowadź opis zdjęcia tutaj

#include <stdio.h>
#include <math.h>
#include <string.h>


#define SAMPLE_LENGTH 1000

float stddev(float data[], int len);
float mean(float data[], int len);
void thresholding(float y[], int signals[], int lag, float threshold, float influence);


void thresholding(float y[], int signals[], int lag, float threshold, float influence) {
    memset(signals, 0, sizeof(float) * SAMPLE_LENGTH);
    float filteredY[SAMPLE_LENGTH];
    memcpy(filteredY, y, sizeof(float) * SAMPLE_LENGTH);
    float avgFilter[SAMPLE_LENGTH];
    float stdFilter[SAMPLE_LENGTH];

    avgFilter[lag - 1] = mean(y, lag);
    stdFilter[lag - 1] = stddev(y, lag);

    for (int i = lag; i < SAMPLE_LENGTH; i++) {
        if (fabsf(y[i] - avgFilter[i-1]) > threshold * stdFilter[i-1]) {
            if (y[i] > avgFilter[i-1]) {
                signals[i] = 1;
            } else {
                signals[i] = -1;
            }
            filteredY[i] = influence * y[i] + (1 - influence) * filteredY[i-1];
        } else {
            signals[i] = 0;
        }
        avgFilter[i] = mean(filteredY + i-lag, lag);
        stdFilter[i] = stddev(filteredY + i-lag, lag);
    }
}

float mean(float data[], int len) {
    float sum = 0.0, mean = 0.0;

    int i;
    for(i=0; i<len; ++i) {
        sum += data[i];
    }

    mean = sum/len;
    return mean;


}

float stddev(float data[], int len) {
    float the_mean = mean(data, len);
    float standardDeviation = 0.0;

    int i;
    for(i=0; i<len; ++i) {
        standardDeviation += pow(data[i] - the_mean, 2);
    }

    return sqrt(standardDeviation/len);
}

int main() {
    printf("Hello, World!\n");
    int lag = 100;
    float threshold = 5;
    float influence = 0;
    float y[]=  {1,1,1.1,1,0.9,1,1,1.1,1,0.9,1,1.1,1,1,0.9,1,1,1.1,1,1,1,1,1.1,0.9,1,1.1,1,1,0.9,
  ....
1,1.1,1,1,1.1,1,0.8,0.9,1,1.2,0.9,1,1,1.1,1.2,1,1.5,1,3,2,5,3,2,1,1,1,0.9,1,1,3,       2.6,4,3,3.2,2,1,1,0.8,4,4,2,2.5,1,1,1,1.2,1,1.5,1,3,2,5,3,2,1,1,1,0.9,1,1,3,
       2.6,4,3,3.2,2,1,1,0.8,4,4,2,2.5,1,1,1}

    int signal[SAMPLE_LENGTH];

    thresholding(y, signal,  lag, threshold, influence);

    return 0;
}

Jej wynik to wpływ = 0

wprowadź opis zdjęcia tutaj

Nie wspaniale, ale tutaj z wpływem = 1

wprowadź opis zdjęcia tutaj

co jest bardzo dobre.

DavidC
źródło
5

Oto rzeczywista implementacja Java oparta na odpowiedzi Groovy opublikowanej wcześniej . (Wiem, że są już publikowane implementacje Groovy i Kotlin, ale dla kogoś takiego jak ja, który tylko napisał Javę, trudno jest znaleźć sposób na konwersję między innymi językami i Javą).

(Wyniki pasują do wykresów innych osób)

Implementacja algorytmu

import java.util.ArrayList;
import java.util.Collections;
import java.util.HashMap;
import java.util.List;

import org.apache.commons.math3.stat.descriptive.SummaryStatistics;

public class SignalDetector {

    public HashMap<String, List> analyzeDataForSignals(List<Double> data, int lag, Double threshold, Double influence) {

        // init stats instance
        SummaryStatistics stats = new SummaryStatistics();

        // the results (peaks, 1 or -1) of our algorithm
        List<Integer> signals = new ArrayList<Integer>(Collections.nCopies(data.size(), 0));

        // filter out the signals (peaks) from our original list (using influence arg)
        List<Double> filteredData = new ArrayList<Double>(data);

        // the current average of the rolling window
        List<Double> avgFilter = new ArrayList<Double>(Collections.nCopies(data.size(), 0.0d));

        // the current standard deviation of the rolling window
        List<Double> stdFilter = new ArrayList<Double>(Collections.nCopies(data.size(), 0.0d));

        // init avgFilter and stdFilter
        for (int i = 0; i < lag; i++) {
            stats.addValue(data.get(i));
        }
        avgFilter.set(lag - 1, stats.getMean());
        stdFilter.set(lag - 1, Math.sqrt(stats.getPopulationVariance())); // getStandardDeviation() uses sample variance
        stats.clear();

        // loop input starting at end of rolling window
        for (int i = lag; i < data.size(); i++) {

            // if the distance between the current value and average is enough standard deviations (threshold) away
            if (Math.abs((data.get(i) - avgFilter.get(i - 1))) > threshold * stdFilter.get(i - 1)) {

                // this is a signal (i.e. peak), determine if it is a positive or negative signal
                if (data.get(i) > avgFilter.get(i - 1)) {
                    signals.set(i, 1);
                } else {
                    signals.set(i, -1);
                }

                // filter this signal out using influence
                filteredData.set(i, (influence * data.get(i)) + ((1 - influence) * filteredData.get(i - 1)));
            } else {
                // ensure this signal remains a zero
                signals.set(i, 0);
                // ensure this value is not filtered
                filteredData.set(i, data.get(i));
            }

            // update rolling average and deviation
            for (int j = i - lag; j < i; j++) {
                stats.addValue(filteredData.get(j));
            }
            avgFilter.set(i, stats.getMean());
            stdFilter.set(i, Math.sqrt(stats.getPopulationVariance()));
            stats.clear();
        }

        HashMap<String, List> returnMap = new HashMap<String, List>();
        returnMap.put("signals", signals);
        returnMap.put("filteredData", filteredData);
        returnMap.put("avgFilter", avgFilter);
        returnMap.put("stdFilter", stdFilter);

        return returnMap;

    } // end
}

Główna metoda

import java.text.DecimalFormat;
import java.util.ArrayList;
import java.util.Arrays;
import java.util.HashMap;
import java.util.List;

public class Main {

    public static void main(String[] args) throws Exception {
        DecimalFormat df = new DecimalFormat("#0.000");

        ArrayList<Double> data = new ArrayList<Double>(Arrays.asList(1d, 1d, 1.1d, 1d, 0.9d, 1d, 1d, 1.1d, 1d, 0.9d, 1d,
                1.1d, 1d, 1d, 0.9d, 1d, 1d, 1.1d, 1d, 1d, 1d, 1d, 1.1d, 0.9d, 1d, 1.1d, 1d, 1d, 0.9d, 1d, 1.1d, 1d, 1d,
                1.1d, 1d, 0.8d, 0.9d, 1d, 1.2d, 0.9d, 1d, 1d, 1.1d, 1.2d, 1d, 1.5d, 1d, 3d, 2d, 5d, 3d, 2d, 1d, 1d, 1d,
                0.9d, 1d, 1d, 3d, 2.6d, 4d, 3d, 3.2d, 2d, 1d, 1d, 0.8d, 4d, 4d, 2d, 2.5d, 1d, 1d, 1d));

        SignalDetector signalDetector = new SignalDetector();
        int lag = 30;
        double threshold = 5;
        double influence = 0;

        HashMap<String, List> resultsMap = signalDetector.analyzeDataForSignals(data, lag, threshold, influence);
        // print algorithm params
        System.out.println("lag: " + lag + "\t\tthreshold: " + threshold + "\t\tinfluence: " + influence);

        System.out.println("Data size: " + data.size());
        System.out.println("Signals size: " + resultsMap.get("signals").size());

        // print data
        System.out.print("Data:\t\t");
        for (double d : data) {
            System.out.print(df.format(d) + "\t");
        }
        System.out.println();

        // print signals
        System.out.print("Signals:\t");
        List<Integer> signalsList = resultsMap.get("signals");
        for (int i : signalsList) {
            System.out.print(df.format(i) + "\t");
        }
        System.out.println();

        // print filtered data
        System.out.print("Filtered Data:\t");
        List<Double> filteredDataList = resultsMap.get("filteredData");
        for (double d : filteredDataList) {
            System.out.print(df.format(d) + "\t");
        }
        System.out.println();

        // print running average
        System.out.print("Avg Filter:\t");
        List<Double> avgFilterList = resultsMap.get("avgFilter");
        for (double d : avgFilterList) {
            System.out.print(df.format(d) + "\t");
        }
        System.out.println();

        // print running std
        System.out.print("Std filter:\t");
        List<Double> stdFilterList = resultsMap.get("stdFilter");
        for (double d : stdFilterList) {
            System.out.print(df.format(d) + "\t");
        }
        System.out.println();

        System.out.println();
        for (int i = 0; i < signalsList.size(); i++) {
            if (signalsList.get(i) != 0) {
                System.out.println("Point " + i + " gave signal " + signalsList.get(i));
            }
        }
    }
}

Wyniki

lag: 30     threshold: 5.0      influence: 0.0
Data size: 74
Signals size: 74
Data:           1.000   1.000   1.100   1.000   0.900   1.000   1.000   1.100   1.000   0.900   1.000   1.100   1.000   1.000   0.900   1.000   1.000   1.100   1.000   1.000   1.000   1.000   1.100   0.900   1.000   1.100   1.000   1.000   0.900   1.000   1.100   1.000   1.000   1.100   1.000   0.800   0.900   1.000   1.200   0.900   1.000   1.000   1.100   1.200   1.000   1.500   1.000   3.000   2.000   5.000   3.000   2.000   1.000   1.000   1.000   0.900   1.000   1.000   3.000   2.600   4.000   3.000   3.200   2.000   1.000   1.000   0.800   4.000   4.000   2.000   2.500   1.000   1.000   1.000   
Signals:        0.000   0.000   0.000   0.000   0.000   0.000   0.000   0.000   0.000   0.000   0.000   0.000   0.000   0.000   0.000   0.000   0.000   0.000   0.000   0.000   0.000   0.000   0.000   0.000   0.000   0.000   0.000   0.000   0.000   0.000   0.000   0.000   0.000   0.000   0.000   0.000   0.000   0.000   0.000   0.000   0.000   0.000   0.000   0.000   0.000   1.000   0.000   1.000   1.000   1.000   1.000   1.000   0.000   0.000   0.000   0.000   0.000   0.000   1.000   1.000   1.000   1.000   1.000   1.000   0.000   0.000   0.000   1.000   1.000   1.000   1.000   0.000   0.000   0.000   
Filtered Data:  1.000   1.000   1.100   1.000   0.900   1.000   1.000   1.100   1.000   0.900   1.000   1.100   1.000   1.000   0.900   1.000   1.000   1.100   1.000   1.000   1.000   1.000   1.100   0.900   1.000   1.100   1.000   1.000   0.900   1.000   1.100   1.000   1.000   1.100   1.000   0.800   0.900   1.000   1.200   0.900   1.000   1.000   1.100   1.200   1.000   1.000   1.000   1.000   1.000   1.000   1.000   1.000   1.000   1.000   1.000   0.900   1.000   1.000   1.000   1.000   1.000   1.000   1.000   1.000   1.000   1.000   0.800   0.800   0.800   0.800   0.800   1.000   1.000   1.000   
Avg Filter:     0.000   0.000   0.000   0.000   0.000   0.000   0.000   0.000   0.000   0.000   0.000   0.000   0.000   0.000   0.000   0.000   0.000   0.000   0.000   0.000   0.000   0.000   0.000   0.000   0.000   0.000   0.000   0.000   0.000   1.003   1.003   1.007   1.007   1.003   1.007   1.010   1.003   1.000   0.997   1.003   1.003   1.003   1.000   1.003   1.010   1.013   1.013   1.013   1.010   1.010   1.010   1.010   1.010   1.007   1.010   1.010   1.003   1.003   1.003   1.007   1.007   1.003   1.003   1.003   1.000   1.000   1.007   1.003   0.997   0.983   0.980   0.973   0.973   0.970   
Std filter:     0.000   0.000   0.000   0.000   0.000   0.000   0.000   0.000   0.000   0.000   0.000   0.000   0.000   0.000   0.000   0.000   0.000   0.000   0.000   0.000   0.000   0.000   0.000   0.000   0.000   0.000   0.000   0.000   0.000   0.060   0.060   0.063   0.063   0.060   0.063   0.060   0.071   0.073   0.071   0.080   0.080   0.080   0.077   0.080   0.087   0.085   0.085   0.085   0.083   0.083   0.083   0.083   0.083   0.081   0.079   0.079   0.080   0.080   0.080   0.077   0.077   0.075   0.075   0.075   0.073   0.073   0.063   0.071   0.080   0.078   0.083   0.089   0.089   0.086   

Point 45 gave signal 1
Point 47 gave signal 1
Point 48 gave signal 1
Point 49 gave signal 1
Point 50 gave signal 1
Point 51 gave signal 1
Point 58 gave signal 1
Point 59 gave signal 1
Point 60 gave signal 1
Point 61 gave signal 1
Point 62 gave signal 1
Point 63 gave signal 1
Point 67 gave signal 1
Point 68 gave signal 1
Point 69 gave signal 1
Point 70 gave signal 1

Wykresy przedstawiające dane i wyniki wykonania Java

takanuva15
źródło
5

Dodatek 1 do oryginalnej odpowiedzi: Matlabi Rtłumaczenia

Kod Matlab

function [signals,avgFilter,stdFilter] = ThresholdingAlgo(y,lag,threshold,influence)
% Initialise signal results
signals = zeros(length(y),1);
% Initialise filtered series
filteredY = y(1:lag+1);
% Initialise filters
avgFilter(lag+1,1) = mean(y(1:lag+1));
stdFilter(lag+1,1) = std(y(1:lag+1));
% Loop over all datapoints y(lag+2),...,y(t)
for i=lag+2:length(y)
    % If new value is a specified number of deviations away
    if abs(y(i)-avgFilter(i-1)) > threshold*stdFilter(i-1)
        if y(i) > avgFilter(i-1)
            % Positive signal
            signals(i) = 1;
        else
            % Negative signal
            signals(i) = -1;
        end
        % Make influence lower
        filteredY(i) = influence*y(i)+(1-influence)*filteredY(i-1);
    else
        % No signal
        signals(i) = 0;
        filteredY(i) = y(i);
    end
    % Adjust the filters
    avgFilter(i) = mean(filteredY(i-lag:i));
    stdFilter(i) = std(filteredY(i-lag:i));
end
% Done, now return results
end

Przykład:

% Data
y = [1 1 1.1 1 0.9 1 1 1.1 1 0.9 1 1.1 1 1 0.9 1 1 1.1 1 1,...
    1 1 1.1 0.9 1 1.1 1 1 0.9 1 1.1 1 1 1.1 1 0.8 0.9 1 1.2 0.9 1,...
    1 1.1 1.2 1 1.5 1 3 2 5 3 2 1 1 1 0.9 1,...
    1 3 2.6 4 3 3.2 2 1 1 0.8 4 4 2 2.5 1 1 1];

% Settings
lag = 30;
threshold = 5;
influence = 0;

% Get results
[signals,avg,dev] = ThresholdingAlgo(y,lag,threshold,influence);

figure; subplot(2,1,1); hold on;
x = 1:length(y); ix = lag+1:length(y);
area(x(ix),avg(ix)+threshold*dev(ix),'FaceColor',[0.9 0.9 0.9],'EdgeColor','none');
area(x(ix),avg(ix)-threshold*dev(ix),'FaceColor',[1 1 1],'EdgeColor','none');
plot(x(ix),avg(ix),'LineWidth',1,'Color','cyan','LineWidth',1.5);
plot(x(ix),avg(ix)+threshold*dev(ix),'LineWidth',1,'Color','green','LineWidth',1.5);
plot(x(ix),avg(ix)-threshold*dev(ix),'LineWidth',1,'Color','green','LineWidth',1.5);
plot(1:length(y),y,'b');
subplot(2,1,2);
stairs(signals,'r','LineWidth',1.5); ylim([-1.5 1.5]);

Kod R.

ThresholdingAlgo <- function(y,lag,threshold,influence) {
  signals <- rep(0,length(y))
  filteredY <- y[0:lag]
  avgFilter <- NULL
  stdFilter <- NULL
  avgFilter[lag] <- mean(y[0:lag], na.rm=TRUE)
  stdFilter[lag] <- sd(y[0:lag], na.rm=TRUE)
  for (i in (lag+1):length(y)){
    if (abs(y[i]-avgFilter[i-1]) > threshold*stdFilter[i-1]) {
      if (y[i] > avgFilter[i-1]) {
        signals[i] <- 1;
      } else {
        signals[i] <- -1;
      }
      filteredY[i] <- influence*y[i]+(1-influence)*filteredY[i-1]
    } else {
      signals[i] <- 0
      filteredY[i] <- y[i]
    }
    avgFilter[i] <- mean(filteredY[(i-lag):i], na.rm=TRUE)
    stdFilter[i] <- sd(filteredY[(i-lag):i], na.rm=TRUE)
  }
  return(list("signals"=signals,"avgFilter"=avgFilter,"stdFilter"=stdFilter))
}

Przykład:

# Data
y <- c(1,1,1.1,1,0.9,1,1,1.1,1,0.9,1,1.1,1,1,0.9,1,1,1.1,1,1,1,1,1.1,0.9,1,1.1,1,1,0.9,
       1,1.1,1,1,1.1,1,0.8,0.9,1,1.2,0.9,1,1,1.1,1.2,1,1.5,1,3,2,5,3,2,1,1,1,0.9,1,1,3,
       2.6,4,3,3.2,2,1,1,0.8,4,4,2,2.5,1,1,1)

lag       <- 30
threshold <- 5
influence <- 0

# Run algo with lag = 30, threshold = 5, influence = 0
result <- ThresholdingAlgo(y,lag,threshold,influence)

# Plot result
par(mfrow = c(2,1),oma = c(2,2,0,0) + 0.1,mar = c(0,0,2,1) + 0.2)
plot(1:length(y),y,type="l",ylab="",xlab="") 
lines(1:length(y),result$avgFilter,type="l",col="cyan",lwd=2)
lines(1:length(y),result$avgFilter+threshold*result$stdFilter,type="l",col="green",lwd=2)
lines(1:length(y),result$avgFilter-threshold*result$stdFilter,type="l",col="green",lwd=2)
plot(result$signals,type="S",col="red",ylab="",xlab="",ylim=c(-1.5,1.5),lwd=2)

Ten kod (oba języki) da następujący wynik dla danych pierwotnego pytania:

Przykład progowania z kodu Matlab


Dodatek 2 do oryginalnej odpowiedzi: Matlabkod demonstracyjny

(kliknij, aby utworzyć dane)

Demo Matlaba

function [] = RobustThresholdingDemo()

%% SPECIFICATIONS
lag         = 5;       % lag for the smoothing
threshold   = 3.5;     % number of st.dev. away from the mean to signal
influence   = 0.3;     % when signal: how much influence for new data? (between 0 and 1)
                       % 1 is normal influence, 0.5 is half      
%% START DEMO
DemoScreen(30,lag,threshold,influence);

end

function [signals,avgFilter,stdFilter] = ThresholdingAlgo(y,lag,threshold,influence)
signals = zeros(length(y),1);
filteredY = y(1:lag+1);
avgFilter(lag+1,1) = mean(y(1:lag+1));
stdFilter(lag+1,1) = std(y(1:lag+1));
for i=lag+2:length(y)
    if abs(y(i)-avgFilter(i-1)) > threshold*stdFilter(i-1)
        if y(i) > avgFilter(i-1)
            signals(i) = 1;
        else
            signals(i) = -1;
        end
        filteredY(i) = influence*y(i)+(1-influence)*filteredY(i-1);
    else
        signals(i) = 0;
        filteredY(i) = y(i);
    end
    avgFilter(i) = mean(filteredY(i-lag:i));
    stdFilter(i) = std(filteredY(i-lag:i));
end
end

% Demo screen function
function [] = DemoScreen(n,lag,threshold,influence)
figure('Position',[200 100,1000,500]);
subplot(2,1,1);
title(sprintf(['Draw data points (%.0f max)      [settings: lag = %.0f, '...
    'threshold = %.2f, influence = %.2f]'],n,lag,threshold,influence));
ylim([0 5]); xlim([0 50]);
H = gca; subplot(2,1,1);
set(H, 'YLimMode', 'manual'); set(H, 'XLimMode', 'manual');
set(H, 'YLim', get(H,'YLim')); set(H, 'XLim', get(H,'XLim'));
xg = []; yg = [];
for i=1:n
    try
        [xi,yi] = ginput(1);
    catch
        return;
    end
    xg = [xg xi]; yg = [yg yi];
    if i == 1
        subplot(2,1,1); hold on;
        plot(H, xg(i),yg(i),'r.'); 
        text(xg(i),yg(i),num2str(i),'FontSize',7);
    end
    if length(xg) > lag
        [signals,avg,dev] = ...
            ThresholdingAlgo(yg,lag,threshold,influence);
        area(xg(lag+1:end),avg(lag+1:end)+threshold*dev(lag+1:end),...
            'FaceColor',[0.9 0.9 0.9],'EdgeColor','none');
        area(xg(lag+1:end),avg(lag+1:end)-threshold*dev(lag+1:end),...
            'FaceColor',[1 1 1],'EdgeColor','none');
        plot(xg(lag+1:end),avg(lag+1:end),'LineWidth',1,'Color','cyan');
        plot(xg(lag+1:end),avg(lag+1:end)+threshold*dev(lag+1:end),...
            'LineWidth',1,'Color','green');
        plot(xg(lag+1:end),avg(lag+1:end)-threshold*dev(lag+1:end),...
            'LineWidth',1,'Color','green');
        subplot(2,1,2); hold on; title('Signal output');
        stairs(xg(lag+1:end),signals(lag+1:end),'LineWidth',2,'Color','blue');
        ylim([-2 2]); xlim([0 50]); hold off;
    end
    subplot(2,1,1); hold on;
    for j=2:i
        plot(xg([j-1:j]),yg([j-1:j]),'r'); plot(H,xg(j),yg(j),'r.');
        text(xg(j),yg(j),num2str(j),'FontSize',7);
    end
end
end

Jean-Paul
źródło
4

Oto moja próba stworzenia rozwiązania Ruby dla „Wygładzonego z-score algo” z zaakceptowanej odpowiedzi:

module ThresholdingAlgoMixin
  def mean(array)
    array.reduce(&:+) / array.size.to_f
  end

  def stddev(array)
    array_mean = mean(array)
    Math.sqrt(array.reduce(0.0) { |a, b| a.to_f + ((b.to_f - array_mean) ** 2) } / array.size.to_f)
  end

  def thresholding_algo(lag: 5, threshold: 3.5, influence: 0.5)
    return nil if size < lag * 2
    Array.new(size, 0).tap do |signals|
      filtered = Array.new(self)

      initial_slice = take(lag)
      avg_filter = Array.new(lag - 1, 0.0) + [mean(initial_slice)]
      std_filter = Array.new(lag - 1, 0.0) + [stddev(initial_slice)]
      (lag..size-1).each do |idx|
        prev = idx - 1
        if (fetch(idx) - avg_filter[prev]).abs > threshold * std_filter[prev]
          signals[idx] = fetch(idx) > avg_filter[prev] ? 1 : -1
          filtered[idx] = (influence * fetch(idx)) + ((1-influence) * filtered[prev])
        end

        filtered_slice = filtered[idx-lag..prev]
        avg_filter[idx] = mean(filtered_slice)
        std_filter[idx] = stddev(filtered_slice)
      end
    end
  end
end

I przykładowe użycie:

test_data = [
  1, 1, 1.1, 1, 0.9, 1, 1, 1.1, 1, 0.9, 1, 1.1, 1, 1, 0.9, 1,
  1, 1.1, 1, 1, 1, 1, 1.1, 0.9, 1, 1.1, 1, 1, 0.9, 1, 1.1, 1,
  1, 1.1, 1, 0.8, 0.9, 1, 1.2, 0.9, 1, 1, 1.1, 1.2, 1, 1.5,
  1, 3, 2, 5, 3, 2, 1, 1, 1, 0.9, 1, 1, 3, 2.6, 4, 3, 3.2, 2,
  1, 1, 0.8, 4, 4, 2, 2.5, 1, 1, 1
].extend(ThresholdingAlgoMixin)

puts test_data.thresholding_algo.inspect

# Output: [
#   0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
#   0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, -1, 0, 0, 0,
#   0, 0, 0, 0, 0, 0, 1, 0, 1, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 1,
#   1, 1, 0, 0, 0, -1, -1, 0, 0, 0, 0, 0, 0, 0, 0
# ]
Kimmo Lehto
źródło
Wspaniale, dziękuję za udostępnienie! Dodam cię do listy. Upewnij się, że dla aplikacji czasu rzeczywistego utworzysz osobną funkcję aktualizującą sygnały, gdy nadejdzie nowy punkt danych (zamiast zapętlać wszystkie punkty danych za każdym razem).
Jean-Paul,
4

Wersja iteracyjna w python / numpy dla odpowiedzi https://stackoverflow.com/a/22640362/6029703 jest tutaj. Ten kod jest szybszy niż średnia obliczeniowa i odchylenie standardowe przy każdym opóźnieniu dla dużych danych (100000+).

def peak_detection_smoothed_zscore_v2(x, lag, threshold, influence):
    '''
    iterative smoothed z-score algorithm
    Implementation of algorithm from https://stackoverflow.com/a/22640362/6029703
    '''
    import numpy as np
    labels = np.zeros(len(x))
    filtered_y = np.array(x)
    avg_filter = np.zeros(len(x))
    std_filter = np.zeros(len(x))
    var_filter = np.zeros(len(x))

    avg_filter[lag - 1] = np.mean(x[0:lag])
    std_filter[lag - 1] = np.std(x[0:lag])
    var_filter[lag - 1] = np.var(x[0:lag])
    for i in range(lag, len(x)):
        if abs(x[i] - avg_filter[i - 1]) > threshold * std_filter[i - 1]:
            if x[i] > avg_filter[i - 1]:
                labels[i] = 1
            else:
                labels[i] = -1
            filtered_y[i] = influence * x[i] + (1 - influence) * filtered_y[i - 1]
        else:
            labels[i] = 0
            filtered_y[i] = x[i]
        # update avg, var, std
        avg_filter[i] = avg_filter[i - 1] + 1. / lag * (filtered_y[i] - filtered_y[i - lag])
        var_filter[i] = var_filter[i - 1] + 1. / lag * ((filtered_y[i] - avg_filter[i - 1]) ** 2 - (
            filtered_y[i - lag] - avg_filter[i - 1]) ** 2 - (filtered_y[i] - filtered_y[i - lag]) ** 2 / lag)
        std_filter[i] = np.sqrt(var_filter[i])

    return dict(signals=labels,
                avgFilter=avg_filter,
                stdFilter=std_filter)
Tranfer Will
źródło
4

Pomyślałem, że zapewnię moją implementację algorytmu Julii dla innych. Istotę można znaleźć tutaj

using Statistics
using Plots
function SmoothedZscoreAlgo(y, lag, threshold, influence)
    # Julia implimentation of http://stackoverflow.com/a/22640362/6029703
    n = length(y)
    signals = zeros(n) # init signal results
    filteredY = copy(y) # init filtered series
    avgFilter = zeros(n) # init average filter
    stdFilter = zeros(n) # init std filter
    avgFilter[lag - 1] = mean(y[1:lag]) # init first value
    stdFilter[lag - 1] = std(y[1:lag]) # init first value

    for i in range(lag, stop=n-1)
        if abs(y[i] - avgFilter[i-1]) > threshold*stdFilter[i-1]
            if y[i] > avgFilter[i-1]
                signals[i] += 1 # postive signal
            else
                signals[i] += -1 # negative signal
            end
            # Make influence lower
            filteredY[i] = influence*y[i] + (1-influence)*filteredY[i-1]
        else
            signals[i] = 0
            filteredY[i] = y[i]
        end
        avgFilter[i] = mean(filteredY[i-lag+1:i])
        stdFilter[i] = std(filteredY[i-lag+1:i])
    end
    return (signals = signals, avgFilter = avgFilter, stdFilter = stdFilter)
end


# Data
y = [1,1,1.1,1,0.9,1,1,1.1,1,0.9,1,1.1,1,1,0.9,1,1,1.1,1,1,1,1,1.1,0.9,1,1.1,1,1,0.9,
       1,1.1,1,1,1.1,1,0.8,0.9,1,1.2,0.9,1,1,1.1,1.2,1,1.5,1,3,2,5,3,2,1,1,1,0.9,1,1,3,
       2.6,4,3,3.2,2,1,1,0.8,4,4,2,2.5,1,1,1]

# Settings: lag = 30, threshold = 5, influence = 0
lag = 30
threshold = 5
influence = 0

results = SmoothedZscoreAlgo(y, lag, threshold, influence)
upper_bound = results[:avgFilter] + threshold * results[:stdFilter]
lower_bound = results[:avgFilter] - threshold * results[:stdFilter]
x = 1:length(y)

yplot = plot(x,y,color="blue", label="Y",legend=:topleft)
yplot = plot!(x,upper_bound, color="green", label="Upper Bound",legend=:topleft)
yplot = plot!(x,results[:avgFilter], color="cyan", label="Average Filter",legend=:topleft)
yplot = plot!(x,lower_bound, color="green", label="Lower Bound",legend=:topleft)
signalplot = plot(x,results[:signals],color="red",label="Signals",legend=:topleft)
plot(yplot,signalplot,layout=(2,1),legend=:topleft)

Wyniki

Matt Camp
źródło
3

Oto implementacja wygładzonego algorytmu Z-score Groovy (Java) ( patrz odpowiedź powyżej ).

/**
 * "Smoothed zero-score alogrithm" shamelessly copied from https://stackoverflow.com/a/22640362/6029703
 *  Uses a rolling mean and a rolling deviation (separate) to identify peaks in a vector
 *
 * @param y - The input vector to analyze
 * @param lag - The lag of the moving window (i.e. how big the window is)
 * @param threshold - The z-score at which the algorithm signals (i.e. how many standard deviations away from the moving mean a peak (or signal) is)
 * @param influence - The influence (between 0 and 1) of new signals on the mean and standard deviation (how much a peak (or signal) should affect other values near it)
 * @return - The calculated averages (avgFilter) and deviations (stdFilter), and the signals (signals)
 */

public HashMap<String, List<Object>> thresholdingAlgo(List<Double> y, Long lag, Double threshold, Double influence) {
    //init stats instance
    SummaryStatistics stats = new SummaryStatistics()

    //the results (peaks, 1 or -1) of our algorithm
    List<Integer> signals = new ArrayList<Integer>(Collections.nCopies(y.size(), 0))
    //filter out the signals (peaks) from our original list (using influence arg)
    List<Double> filteredY = new ArrayList<Double>(y)
    //the current average of the rolling window
    List<Double> avgFilter = new ArrayList<Double>(Collections.nCopies(y.size(), 0.0d))
    //the current standard deviation of the rolling window
    List<Double> stdFilter = new ArrayList<Double>(Collections.nCopies(y.size(), 0.0d))
    //init avgFilter and stdFilter
    (0..lag-1).each { stats.addValue(y[it as int]) }
    avgFilter[lag - 1 as int] = stats.getMean()
    stdFilter[lag - 1 as int] = Math.sqrt(stats.getPopulationVariance()) //getStandardDeviation() uses sample variance (not what we want)
    stats.clear()
    //loop input starting at end of rolling window
    (lag..y.size()-1).each { i ->
        //if the distance between the current value and average is enough standard deviations (threshold) away
        if (Math.abs((y[i as int] - avgFilter[i - 1 as int]) as Double) > threshold * stdFilter[i - 1 as int]) {
            //this is a signal (i.e. peak), determine if it is a positive or negative signal
            signals[i as int] = (y[i as int] > avgFilter[i - 1 as int]) ? 1 : -1
            //filter this signal out using influence
            filteredY[i as int] = (influence * y[i as int]) + ((1-influence) * filteredY[i - 1 as int])
        } else {
            //ensure this signal remains a zero
            signals[i as int] = 0
            //ensure this value is not filtered
            filteredY[i as int] = y[i as int]
        }
        //update rolling average and deviation
        (i - lag..i-1).each { stats.addValue(filteredY[it as int] as Double) }
        avgFilter[i as int] = stats.getMean()
        stdFilter[i as int] = Math.sqrt(stats.getPopulationVariance()) //getStandardDeviation() uses sample variance (not what we want)
        stats.clear()
    }

    return [
        signals  : signals,
        avgFilter: avgFilter,
        stdFilter: stdFilter
    ]
}

Poniżej znajduje się test na tym samym zestawie danych, który daje takie same wyniki jak powyższa implementacja Python / numpy .

    // Data
    def y = [1d, 1d, 1.1d, 1d, 0.9d, 1d, 1d, 1.1d, 1d, 0.9d, 1d, 1.1d, 1d, 1d, 0.9d, 1d, 1d, 1.1d, 1d, 1d,
         1d, 1d, 1.1d, 0.9d, 1d, 1.1d, 1d, 1d, 0.9d, 1d, 1.1d, 1d, 1d, 1.1d, 1d, 0.8d, 0.9d, 1d, 1.2d, 0.9d, 1d,
         1d, 1.1d, 1.2d, 1d, 1.5d, 1d, 3d, 2d, 5d, 3d, 2d, 1d, 1d, 1d, 0.9d, 1d,
         1d, 3d, 2.6d, 4d, 3d, 3.2d, 2d, 1d, 1d, 0.8d, 4d, 4d, 2d, 2.5d, 1d, 1d, 1d]

    // Settings
    def lag = 30
    def threshold = 5
    def influence = 0


    def thresholdingResults = thresholdingAlgo((List<Double>) y, (Long) lag, (Double) threshold, (Double) influence)

    println y.size()
    println thresholdingResults.signals.size()
    println thresholdingResults.signals

    thresholdingResults.signals.eachWithIndex { x, idx ->
        if (x) {
            println y[idx]
        }
    }
JoshuaCWebDeveloper
źródło
3

Oto (nieidiomatyczna) wersja Scala wygładzonego algorytmu Z-score :

/**
  * Smoothed zero-score alogrithm shamelessly copied from https://stackoverflow.com/a/22640362/6029703
  * Uses a rolling mean and a rolling deviation (separate) to identify peaks in a vector
  *
  * @param y - The input vector to analyze
  * @param lag - The lag of the moving window (i.e. how big the window is)
  * @param threshold - The z-score at which the algorithm signals (i.e. how many standard deviations away from the moving mean a peak (or signal) is)
  * @param influence - The influence (between 0 and 1) of new signals on the mean and standard deviation (how much a peak (or signal) should affect other values near it)
  * @return - The calculated averages (avgFilter) and deviations (stdFilter), and the signals (signals)
  */
private def smoothedZScore(y: Seq[Double], lag: Int, threshold: Double, influence: Double): Seq[Int] = {
  val stats = new SummaryStatistics()

  // the results (peaks, 1 or -1) of our algorithm
  val signals = mutable.ArrayBuffer.fill(y.length)(0)

  // filter out the signals (peaks) from our original list (using influence arg)
  val filteredY = y.to[mutable.ArrayBuffer]

  // the current average of the rolling window
  val avgFilter = mutable.ArrayBuffer.fill(y.length)(0d)

  // the current standard deviation of the rolling window
  val stdFilter = mutable.ArrayBuffer.fill(y.length)(0d)

  // init avgFilter and stdFilter
  y.take(lag).foreach(s => stats.addValue(s))

  avgFilter(lag - 1) = stats.getMean
  stdFilter(lag - 1) = Math.sqrt(stats.getPopulationVariance) // getStandardDeviation() uses sample variance (not what we want)

  // loop input starting at end of rolling window
  y.zipWithIndex.slice(lag, y.length - 1).foreach {
    case (s: Double, i: Int) =>
      // if the distance between the current value and average is enough standard deviations (threshold) away
      if (Math.abs(s - avgFilter(i - 1)) > threshold * stdFilter(i - 1)) {
        // this is a signal (i.e. peak), determine if it is a positive or negative signal
        signals(i) = if (s > avgFilter(i - 1)) 1 else -1
        // filter this signal out using influence
        filteredY(i) = (influence * s) + ((1 - influence) * filteredY(i - 1))
      } else {
        // ensure this signal remains a zero
        signals(i) = 0
        // ensure this value is not filtered
        filteredY(i) = s
      }

      // update rolling average and deviation
      stats.clear()
      filteredY.slice(i - lag, i).foreach(s => stats.addValue(s))
      avgFilter(i) = stats.getMean
      stdFilter(i) = Math.sqrt(stats.getPopulationVariance) // getStandardDeviation() uses sample variance (not what we want)
  }

  println(y.length)
  println(signals.length)
  println(signals)

  signals.zipWithIndex.foreach {
    case(x: Int, idx: Int) =>
      if (x == 1) {
        println(idx + " " + y(idx))
      }
  }

  val data =
    y.zipWithIndex.map { case (s: Double, i: Int) => Map("x" -> i, "y" -> s, "name" -> "y", "row" -> "data") } ++
    avgFilter.zipWithIndex.map { case (s: Double, i: Int) => Map("x" -> i, "y" -> s, "name" -> "avgFilter", "row" -> "data") } ++
    avgFilter.zipWithIndex.map { case (s: Double, i: Int) => Map("x" -> i, "y" -> (s - threshold * stdFilter(i)), "name" -> "lower", "row" -> "data") } ++
    avgFilter.zipWithIndex.map { case (s: Double, i: Int) => Map("x" -> i, "y" -> (s + threshold * stdFilter(i)), "name" -> "upper", "row" -> "data") } ++
    signals.zipWithIndex.map { case (s: Int, i: Int) => Map("x" -> i, "y" -> s, "name" -> "signal", "row" -> "signal") }

  Vegas("Smoothed Z")
    .withData(data)
    .mark(Line)
    .encodeX("x", Quant)
    .encodeY("y", Quant)
    .encodeColor(
      field="name",
      dataType=Nominal
    )
    .encodeRow("row", Ordinal)
    .show

  return signals
}

Oto test, który zwraca takie same wyniki jak wersje Python i Groovy:

val y = List(1d, 1d, 1.1d, 1d, 0.9d, 1d, 1d, 1.1d, 1d, 0.9d, 1d, 1.1d, 1d, 1d, 0.9d, 1d, 1d, 1.1d, 1d, 1d,
  1d, 1d, 1.1d, 0.9d, 1d, 1.1d, 1d, 1d, 0.9d, 1d, 1.1d, 1d, 1d, 1.1d, 1d, 0.8d, 0.9d, 1d, 1.2d, 0.9d, 1d,
  1d, 1.1d, 1.2d, 1d, 1.5d, 1d, 3d, 2d, 5d, 3d, 2d, 1d, 1d, 1d, 0.9d, 1d,
  1d, 3d, 2.6d, 4d, 3d, 3.2d, 2d, 1d, 1d, 0.8d, 4d, 4d, 2d, 2.5d, 1d, 1d, 1d)

val lag = 30
val threshold = 5d
val influence = 0d

smoothedZScore(y, lag, threshold, influence)

Tabela wyników Vegas

Gist tutaj

Mike Roberts
źródło
1 oznacza szczyty, -1 oznacza doliny.
Mike Roberts,
3

Potrzebowałem czegoś takiego w moim projekcie na Androida. Pomyślałem, że mogę oddać implementację Kotlina .

/**
* Smoothed zero-score alogrithm shamelessly copied from https://stackoverflow.com/a/22640362/6029703
* Uses a rolling mean and a rolling deviation (separate) to identify peaks in a vector
*
* @param y - The input vector to analyze
* @param lag - The lag of the moving window (i.e. how big the window is)
* @param threshold - The z-score at which the algorithm signals (i.e. how many standard deviations away from the moving mean a peak (or signal) is)
* @param influence - The influence (between 0 and 1) of new signals on the mean and standard deviation (how much a peak (or signal) should affect other values near it)
* @return - The calculated averages (avgFilter) and deviations (stdFilter), and the signals (signals)
*/
fun smoothedZScore(y: List<Double>, lag: Int, threshold: Double, influence: Double): Triple<List<Int>, List<Double>, List<Double>> {
    val stats = SummaryStatistics()
    // the results (peaks, 1 or -1) of our algorithm
    val signals = MutableList<Int>(y.size, { 0 })
    // filter out the signals (peaks) from our original list (using influence arg)
    val filteredY = ArrayList<Double>(y)
    // the current average of the rolling window
    val avgFilter = MutableList<Double>(y.size, { 0.0 })
    // the current standard deviation of the rolling window
    val stdFilter = MutableList<Double>(y.size, { 0.0 })
    // init avgFilter and stdFilter
    y.take(lag).forEach { s -> stats.addValue(s) }
    avgFilter[lag - 1] = stats.mean
    stdFilter[lag - 1] = Math.sqrt(stats.populationVariance) // getStandardDeviation() uses sample variance (not what we want)
    stats.clear()
    //loop input starting at end of rolling window
    (lag..y.size - 1).forEach { i ->
        //if the distance between the current value and average is enough standard deviations (threshold) away
        if (Math.abs(y[i] - avgFilter[i - 1]) > threshold * stdFilter[i - 1]) {
            //this is a signal (i.e. peak), determine if it is a positive or negative signal
            signals[i] = if (y[i] > avgFilter[i - 1]) 1 else -1
            //filter this signal out using influence
            filteredY[i] = (influence * y[i]) + ((1 - influence) * filteredY[i - 1])
        } else {
            //ensure this signal remains a zero
            signals[i] = 0
            //ensure this value is not filtered
            filteredY[i] = y[i]
        }
        //update rolling average and deviation
        (i - lag..i - 1).forEach { stats.addValue(filteredY[it]) }
        avgFilter[i] = stats.getMean()
        stdFilter[i] = Math.sqrt(stats.getPopulationVariance()) //getStandardDeviation() uses sample variance (not what we want)
        stats.clear()
    }
    return Triple(signals, avgFilter, stdFilter)
}

przykładowy projekt z wykresami weryfikacyjnymi można znaleźć na github .

wprowadź opis zdjęcia tutaj

leonardkraemer
źródło
Niesamowite! Dzięki za udostępnienie. W przypadku aplikacji w czasie rzeczywistym należy utworzyć osobną funkcję, która oblicza nowy sygnał dla każdego przychodzącego punktu danych. Nie zapętlaj pełnych danych za każdym razem, gdy nadejdzie nowy punkt danych, byłoby to bardzo nieefektywne :)
Jean-Paul,
1
Dobra uwaga, nie myślałem o tym, ponieważ okna, których używam, nie nakładają się.
leonardkraemer
3

Oto zmieniona wersja algorytmu z-score Fortrana . Jest on zmieniany specjalnie dla detekcji piku (rezonansu) w funkcjach przenoszenia w przestrzeni częstotliwości (każda zmiana ma mały komentarz w kodzie).

Pierwsza modyfikacja ostrzega użytkownika, jeśli w pobliżu dolnej granicy wektora wejściowego występuje rezonans, wskazany przez odchylenie standardowe wyższe niż określony próg (w tym przypadku 10%). Oznacza to po prostu, że sygnał nie jest wystarczająco płaski, aby detekcja prawidłowo zainicjowała filtry.

Druga modyfikacja polega na tym, że tylko najwyższa wartość piku jest dodawana do znalezionych pików. Osiąga się to poprzez porównanie każdej znalezionej wartości szczytowej z wielkością jej (opóźnionych) poprzedników i jej (opóźnionych) następców.

Trzecią zmianą jest uwzględnienie, że piki rezonansowe zwykle wykazują pewną formę symetrii wokół częstotliwości rezonansowej. Zatem naturalne jest obliczanie średniej i standardowej symetrycznie wokół bieżącego punktu danych (a nie tylko dla poprzedników). Powoduje to lepsze zachowanie w zakresie wykrywania pików.

Modyfikacje powodują, że cały sygnał musi być wcześniej znany funkcji, co jest typowym przypadkiem wykrywania rezonansu (coś w rodzaju Matlaba, przykład Jean-Paula, w którym punkty danych są generowane w locie, nie zadziała).

function PeakDetect(y,lag,threshold, influence)
    implicit none
    ! Declaring part
    real, dimension(:), intent(in) :: y
    integer, dimension(size(y)) :: PeakDetect
    real, dimension(size(y)) :: filteredY, avgFilter, stdFilter
    integer :: lag, ii
    real :: threshold, influence

    ! Executing part
    PeakDetect = 0
    filteredY = 0.0
    filteredY(1:lag+1) = y(1:lag+1)
    avgFilter = 0.0
    avgFilter(lag+1) = mean(y(1:2*lag+1))
    stdFilter = 0.0
    stdFilter(lag+1) = std(y(1:2*lag+1))

    if (stdFilter(lag+1)/avgFilter(lag+1)>0.1) then ! If the coefficient of variation exceeds 10%, the signal is too uneven at the start, possibly because of a peak.
        write(unit=*,fmt=1001)
1001        format(1X,'Warning: Peak detection might have failed, as there may be a peak at the edge of the frequency range.',/)
    end if
    do ii = lag+2, size(y)
        if (abs(y(ii) - avgFilter(ii-1)) > threshold * stdFilter(ii-1)) then
            ! Find only the largest outstanding value which is only the one greater than its predecessor and its successor
            if (y(ii) > avgFilter(ii-1) .AND. y(ii) > y(ii-1) .AND. y(ii) > y(ii+1)) then
                PeakDetect(ii) = 1
            end if
            filteredY(ii) = influence * y(ii) + (1 - influence) * filteredY(ii-1)
        else
            filteredY(ii) = y(ii)
        end if
        ! Modified with respect to the original code. Mean and standard deviation are calculted symmetrically around the current point
        avgFilter(ii) = mean(filteredY(ii-lag:ii+lag))
        stdFilter(ii) = std(filteredY(ii-lag:ii+lag))
    end do
end function PeakDetect

real function mean(y)
    !> @brief Calculates the mean of vector y
    implicit none
    ! Declaring part
    real, dimension(:), intent(in) :: y
    integer :: N
    ! Executing part
    N = max(1,size(y))
    mean = sum(y)/N
end function mean

real function std(y)
    !> @brief Calculates the standard deviation of vector y
    implicit none
    ! Declaring part
    real, dimension(:), intent(in) :: y
    integer :: N
    ! Executing part
    N = max(1,size(y))
    std = sqrt((N*dot_product(y,y) - sum(y)**2) / (N*(N-1)))
end function std

W mojej aplikacji algorytm działa jak urok! wprowadź opis zdjęcia tutaj

Chociaż
źródło
3

Jeśli masz dane w tabeli bazy danych, oto wersja SQL prostego algorytmu z-score:

with data_with_zscore as (
    select
        date_time,
        value,
        value / (avg(value) over ()) as pct_of_mean,
        (value - avg(value) over ()) / (stdev(value) over ()) as z_score
    from {{tablename}}  where datetime > '2018-11-26' and datetime < '2018-12-03'
)


-- select all
select * from data_with_zscore 

-- select only points greater than a certain threshold
select * from data_with_zscore where z_score > abs(2)
Ocean Airdrop
źródło
Twój kod robi coś innego niż algorytm, który zaproponowałem. Twoje zapytanie po prostu oblicza wyniki Z ([punkt danych - średnia] / std), ale nie uwzględnia logiki mojego algorytmu, który ignoruje przeszłe sygnały podczas obliczania nowych progów sygnału. Ignorujesz również trzy parametry (opóźnienie, wpływ, próg). Czy możesz zrewidować swoją odpowiedź, aby uwzględnić rzeczywistą logikę?
Jean-Paul
1
Tak masz rację. Na początku myślałem, że uda mi się uciec z powyższą uproszczoną wersją. Od tego czasu wziąłem twoje pełne rozwiązanie i przeniosłem je do C #. Zobacz moją odpowiedź poniżej. Kiedy będę miał więcej czasu, ponownie odwiedzę tę wersję SQL i wprowadzę Twój algorytm. Nawiasem mówiąc, dziękuję za tak świetną odpowiedź i wizualne wyjaśnienie.
Ocean Airdrop,
Nie ma problemu i cieszę się, że algorytm może ci pomóc! Dziękujemy za przesłanie C #, którego wciąż brakowało. Dodam to do listy tłumaczeń!
Jean-Paul
3

Wersja Python, która działa ze strumieniami w czasie rzeczywistym (nie przelicza wszystkich punktów danych po przybyciu każdego nowego punktu danych). Możesz dostosować to, co zwraca funkcja klasy - do moich celów potrzebowałem tylko sygnałów.

import numpy as np

class real_time_peak_detection():
    def __init__(self, array, lag, threshold, influence):
        self.y = list(array)
        self.length = len(self.y)
        self.lag = lag
        self.threshold = threshold
        self.influence = influence
        self.signals = [0] * len(self.y)
        self.filteredY = np.array(self.y).tolist()
        self.avgFilter = [0] * len(self.y)
        self.stdFilter = [0] * len(self.y)
        self.avgFilter[self.lag - 1] = np.mean(self.y[0:self.lag]).tolist()
        self.stdFilter[self.lag - 1] = np.std(self.y[0:self.lag]).tolist()

    def thresholding_algo(self, new_value):
        self.y.append(new_value)
        i = len(self.y) - 1
        self.length = len(self.y)
        if i < self.lag:
            return 0
        elif i == self.lag:
            self.signals = [0] * len(self.y)
            self.filteredY = np.array(self.y).tolist()
            self.avgFilter = [0] * len(self.y)
            self.stdFilter = [0] * len(self.y)
            self.avgFilter[self.lag - 1] = np.mean(self.y[0:self.lag]).tolist()
            self.stdFilter[self.lag - 1] = np.std(self.y[0:self.lag]).tolist()
            return 0

        self.signals += [0]
        self.filteredY += [0]
        self.avgFilter += [0]
        self.stdFilter += [0]

        if abs(self.y[i] - self.avgFilter[i - 1]) > self.threshold * self.stdFilter[i - 1]:
            if self.y[i] > self.avgFilter[i - 1]:
                self.signals[i] = 1
            else:
                self.signals[i] = -1

            self.filteredY[i] = self.influence * self.y[i] + (1 - self.influence) * self.filteredY[i - 1]
            self.avgFilter[i] = np.mean(self.filteredY[(i - self.lag):i])
            self.stdFilter[i] = np.std(self.filteredY[(i - self.lag):i])
        else:
            self.signals[i] = 0
            self.filteredY[i] = self.y[i]
            self.avgFilter[i] = np.mean(self.filteredY[(i - self.lag):i])
            self.stdFilter[i] = np.std(self.filteredY[(i - self.lag):i])

        return self.signals[i]
delica
źródło
Dziękujemy za opublikowanie, dodałem twoje tłumaczenie do listy.
Jean-Paul,
3

Pozwoliłem sobie na stworzenie wersji javascript. Może to być pomocne. Javascript powinien być bezpośrednią transkrypcją Pseudokodu podanego powyżej. Dostępne jako pakiet npm i repozytorium github:

Tłumaczenie Javascript:

// javascript port of: /programming/22583391/peak-signal-detection-in-realtime-timeseries-data/48895639#48895639

function sum(a) {
    return a.reduce((acc, val) => acc + val)
}

function mean(a) {
    return sum(a) / a.length
}

function stddev(arr) {
    const arr_mean = mean(arr)
    const r = function(acc, val) {
        return acc + ((val - arr_mean) * (val - arr_mean))
    }
    return Math.sqrt(arr.reduce(r, 0.0) / arr.length)
}

function smoothed_z_score(y, params) {
    var p = params || {}
    // init cooefficients
    const lag = p.lag || 5
    const threshold = p.threshold || 3.5
    const influence = p.influece || 0.5

    if (y === undefined || y.length < lag + 2) {
        throw ` ## y data array to short(${y.length}) for given lag of ${lag}`
    }
    //console.log(`lag, threshold, influence: ${lag}, ${threshold}, ${influence}`)

    // init variables
    var signals = Array(y.length).fill(0)
    var filteredY = y.slice(0)
    const lead_in = y.slice(0, lag)
    //console.log("1: " + lead_in.toString())

    var avgFilter = []
    avgFilter[lag - 1] = mean(lead_in)
    var stdFilter = []
    stdFilter[lag - 1] = stddev(lead_in)
    //console.log("2: " + stdFilter.toString())

    for (var i = lag; i < y.length; i++) {
        //console.log(`${y[i]}, ${avgFilter[i-1]}, ${threshold}, ${stdFilter[i-1]}`)
        if (Math.abs(y[i] - avgFilter[i - 1]) > (threshold * stdFilter[i - 1])) {
            if (y[i] > avgFilter[i - 1]) {
                signals[i] = +1 // positive signal
            } else {
                signals[i] = -1 // negative signal
            }
            // make influence lower
            filteredY[i] = influence * y[i] + (1 - influence) * filteredY[i - 1]
        } else {
            signals[i] = 0 // no signal
            filteredY[i] = y[i]
        }

        // adjust the filters
        const y_lag = filteredY.slice(i - lag, i)
        avgFilter[i] = mean(y_lag)
        stdFilter[i] = stddev(y_lag)
    }

    return signals
}

module.exports = smoothed_z_score
Dirk Lüsebrink
źródło
Dziękujemy za opublikowanie tłumaczenia. Dodałem twój kod do twojej odpowiedzi, aby ludzie mogli go szybko zobaczyć. Dodam twoje tłumaczenie do listy.
Jean-Paul,
Do tej pory zaimportowałem jakiś inny algorytm do javascript. Tym razem z numerycznego pyhon, który daje mi większą kontrolę i działa lepiej dla mnie. Pakowane również w npm, a więcej informacji na temat algo z uniwersytetu stanowego w Waszyngtonie można znaleźć na ich stronie jupyter. npmjs.com/package/@joe_six/duarte-watanabe-peak-detection
Dirk Lüsebrink
2

Jeśli wartość graniczna lub inne kryteria zależą od przyszłych wartości, wówczas jedynym rozwiązaniem (bez wehikułu czasu lub innej wiedzy o przyszłych wartościach) jest opóźnienie jakiejkolwiek decyzji do momentu uzyskania wystarczających przyszłych wartości. Jeśli chcesz poziom powyżej średniej, która obejmuje na przykład 20 punktów, musisz poczekać, aż będziesz miał co najmniej 19 punktów przed jakąkolwiek szczytową decyzją, w przeciwnym razie następny nowy punkt może całkowicie zejść z progu 19 punktów temu .

Twój aktualny wykres nie ma żadnych szczytów ... chyba że z góry wiesz, że następnym punktem nie jest 1e99, który po przeskalowaniu wymiaru Y wykresu byłby płaski do tego punktu.

hotpaw2
źródło
Jak powiedziałem wcześniej, możemy założyć, że jeśli wystąpi szczyt, jest on tak duży jak szczyty na zdjęciu i znacznie odbiega od „normalnych” wartości.
Jean-Paul
Jeśli wiesz, jak duże będą szczyty z góry, ustaw wcześniej średnią i / lub próg na wartość poniżej tej wartości.
hotpaw2
1
I tego dokładnie nie wiem.
Jean-Paul
1
Po prostu zaprzeczyłeś sobie i napisałeś, że szczyty mają rozmiar na zdjęciu. Albo to wiesz, albo nie.
hotpaw2
2
Próbuję ci to wyjaśnić. Masz pomysł teraz, prawda? „Jak zidentyfikować znacząco duże szczyty”. Możesz podejść do problemu statystycznie lub za pomocą inteligentnego algorytmu. Z .. As large as in the picturemi chodziło: o podobnych sytuacjach, w których istnieją znaczące szczyty i podstawowy poziom hałasu.
Jean-Paul
2

I tu jest implementacja PHP ZSCORE algo:

<?php
$y = array(1,7,1.1,1,0.9,1,1,1.1,1,0.9,1,1.1,1,1,0.9,1,1,1.1,1,1,1,1,1.1,0.9,1,1.1,1,1,0.9,
       1,1.1,1,1,1.1,1,0.8,0.9,1,1.2,0.9,1,1,1.1,1.2,1,1.5,10,3,2,5,3,2,1,1,1,0.9,1,1,3,
       2.6,4,3,3.2,2,1,1,0.8,4,4,2,2.5,1,1,1);

function mean($data, $start, $len) {
    $avg = 0;
    for ($i = $start; $i < $start+ $len; $i ++)
        $avg += $data[$i];
    return $avg / $len;
}

function stddev($data, $start,$len) {
    $mean = mean($data,$start,$len);
    $dev = 0;
    for ($i = $start; $i < $start+$len; $i++) 
        $dev += (($data[$i] - $mean) * ($data[$i] - $mean));
    return sqrt($dev / $len);
}

function zscore($data, $len, $lag= 20, $threshold = 1, $influence = 1) {

    $signals = array();
    $avgFilter = array();
    $stdFilter = array();
    $filteredY = array();
    $avgFilter[$lag - 1] = mean($data, 0, $lag);
    $stdFilter[$lag - 1] = stddev($data, 0, $lag);

    for ($i = 0; $i < $len; $i++) {
        $filteredY[$i] = $data[$i];
        $signals[$i] = 0;
    }


    for ($i=$lag; $i < $len; $i++) {
        if (abs($data[$i] - $avgFilter[$i-1]) > $threshold * $stdFilter[$lag - 1]) {
            if ($data[$i] > $avgFilter[$i-1]) {
                $signals[$i] = 1;
            }
            else {
                $signals[$i] = -1;
            }
            $filteredY[$i] = $influence * $data[$i] + (1 - $influence) * $filteredY[$i-1];
        } 
        else {
            $signals[$i] = 0;
            $filteredY[$i] = $data[$i];
        }

        $avgFilter[$i] = mean($filteredY, $i - $lag, $lag);
        $stdFilter[$i] = stddev($filteredY, $i - $lag, $lag);
    }
    return $signals;
}

$sig = zscore($y, count($y));

print_r($y); echo "<br><br>";
print_r($sig); echo "<br><br>";

for ($i = 0; $i < count($y); $i++) echo $i. " " . $y[$i]. " ". $sig[$i]."<br>";

?>
radhoo
źródło
Dziękujemy za wysłanie, dodałem twoje tłumaczenie do listy.
Jean-Paul,
1
Jedna uwaga: zważywszy, że algorytm ten będzie stosowany głównie próbek danych, proponuję wdrożyć próbki odchylenie standardowe przez podzielenie ($len - 1)zamiast $lenwstddev()
Jean-Paul
1

Zamiast porównywać maksima ze średnią, można również porównać maksima z sąsiadującymi minimami, w których minima są zdefiniowane tylko powyżej progu hałasu. Jeśli lokalne maksimum jest> 3 razy (lub inny współczynnik ufności) albo sąsiadujące minima, to maksima te są pikiem. Określanie pików jest dokładniejsze w przypadku szerszych ruchomych okien. Nawiasem mówiąc, powyższe wykorzystuje obliczenia wyśrodkowane na środku okna, a nie obliczenia na końcu okna (== lag).

Zauważ, że maksima należy postrzegać jako wzrost sygnału przed i spadek po.

nichole
źródło
1

Funkcja scipy.signal.find_peaks, jak sama nazwa wskazuje, jest do tego przydatna. Ale ważne jest, aby dobrze zrozumieć jego parametry width, threshold, distance a przede wszystkimprominence , aby uzyskać dobry ekstrakcji szczytową.

Zgodnie z moimi testami i dokumentacją, koncepcja wyeksponowania jest „użyteczną koncepcją” do utrzymania dobrych pików i odrzucenia hałaśliwych pików.

Jakie jest znaczenie (topograficzne) ? Jest to „minimalna wysokość niezbędna do zejścia ze szczytu na dowolny wyższy teren” , jak widać tutaj:

Chodzi o to:

Im wyższe znaczenie, tym „ważny” szczyt.

mrk
źródło
1

Obiektowa wersja algorytmu z-score wykorzystująca mordern C +++

template<typename T>
class FindPeaks{
private:
    std::vector<T> m_input_signal;                      // stores input vector
    std::vector<T> m_array_peak_positive;               
    std::vector<T> m_array_peak_negative;               

public:
    FindPeaks(const std::vector<T>& t_input_signal): m_input_signal{t_input_signal}{ }

    void estimate(){
        int lag{5};
        T threshold{ 5 };                                                                                       // set a threshold
        T influence{ 0.5 };                                                                                    // value between 0 to 1, 1 is normal influence and 0.5 is half the influence

        std::vector<T> filtered_signal(m_input_signal.size(), 0.0);                                             // placeholdered for smooth signal, initialie with all zeros
        std::vector<int> signal(m_input_signal.size(), 0);                                                          // vector that stores where the negative and positive located
        std::vector<T> avg_filtered(m_input_signal.size(), 0.0);                                                // moving averages
        std::vector<T> std_filtered(m_input_signal.size(), 0.0);                                                // moving standard deviation

        avg_filtered[lag] = findMean(m_input_signal.begin(), m_input_signal.begin() + lag);                         // pass the iteartor to vector
        std_filtered[lag] = findStandardDeviation(m_input_signal.begin(), m_input_signal.begin() + lag);

        for (size_t iLag = lag + 1; iLag < m_input_signal.size(); ++iLag) {                                         // start index frm 
            if (std::abs(m_input_signal[iLag] - avg_filtered[iLag - 1]) > threshold * std_filtered[iLag - 1]) {     // check if value is above threhold             
                if ((m_input_signal[iLag]) > avg_filtered[iLag - 1]) {
                    signal[iLag] = 1;                                                                               // assign positive signal
                }
                else {
                    signal[iLag] = -1;                                                                                  // assign negative signal
                }
                filtered_signal[iLag] = influence * m_input_signal[iLag] + (1 - influence) * filtered_signal[iLag - 1];        // exponential smoothing
            }
            else {
                signal[iLag] = 0;                                                                                         // no signal
                filtered_signal[iLag] = m_input_signal[iLag];
            }

            avg_filtered[iLag] = findMean(filtered_signal.begin() + (iLag - lag), filtered_signal.begin() + iLag);
            std_filtered[iLag] = findStandardDeviation(filtered_signal.begin() + (iLag - lag), filtered_signal.begin() + iLag);

        }

        for (size_t iSignal = 0; iSignal < m_input_signal.size(); ++iSignal) {
            if (signal[iSignal] == 1) {
                m_array_peak_positive.emplace_back(m_input_signal[iSignal]);                                        // store the positive peaks
            }
            else if (signal[iSignal] == -1) {
                m_array_peak_negative.emplace_back(m_input_signal[iSignal]);                                         // store the negative peaks
            }
        }
        printVoltagePeaks(signal, m_input_signal);

    }

    std::pair< std::vector<T>, std::vector<T> > get_peaks()
    {
        return std::make_pair(m_array_peak_negative, m_array_peak_negative);
    }

};


template<typename T1, typename T2 >
void printVoltagePeaks(std::vector<T1>& m_signal, std::vector<T2>& m_input_signal) {
    std::ofstream output_file("./voltage_peak.csv");
    std::ostream_iterator<T2> output_iterator_voltage(output_file, ",");
    std::ostream_iterator<T1> output_iterator_signal(output_file, ",");
    std::copy(m_input_signal.begin(), m_input_signal.end(), output_iterator_voltage);
    output_file << "\n";
    std::copy(m_signal.begin(), m_signal.end(), output_iterator_signal);
}

template<typename iterator_type>
typename std::iterator_traits<iterator_type>::value_type findMean(iterator_type it, iterator_type end)
{
    /* function that receives iterator to*/
    typename std::iterator_traits<iterator_type>::value_type sum{ 0.0 };
    int counter = 0;
    while (it != end) {
        sum += *(it++);
        counter++;
    }
    return sum / counter;
}

template<typename iterator_type>
typename std::iterator_traits<iterator_type>::value_type findStandardDeviation(iterator_type it, iterator_type end)
{
    auto mean = findMean(it, end);
    typename std::iterator_traits<iterator_type>::value_type sum_squared_error{ 0.0 };
    int counter{ 0 };
    while (it != end) {
        sum_squared_error += std::pow((*(it++) - mean), 2);
        counter++;
    }
    auto standard_deviation = std::sqrt(sum_squared_error / (counter - 1));
    return standard_deviation;
}
Spandyie
źródło
2
Niezłe tłumaczenie. Byłoby nieco ładniejszy jeśli obiekt oszczędza również filtered_signal, signal, avg_filteredi std_filteredzmienne jako prywatne i tylko aktualizuje te tablice raz , kiedy nowy Datapoint przybywa (obecnie pętle kod ponad wszystkich punktów danych za każdym razem to się nazywa). Poprawiłoby to wydajność twojego kodu i jeszcze lepiej odpowiada strukturze OOP.
Jean-Paul