Prosty algorytm do wykrywania wartości odstających online ogólnych szeregów czasowych

Pracuję z dużą ilością szeregów czasowych. Te szeregi czasowe są w zasadzie pomiarami sieci przychodzącymi co 10 minut, a niektóre z nich są okresowe (tj. Przepustowość), a inne nie (tj. Wielkość ruchu routingu).

Chciałbym prosty algorytm do przeprowadzania online „wykrywania wartości odstających”. Zasadniczo chcę zachować w pamięci (lub na dysku) całe dane historyczne dla każdej serii czasowej i chcę wykryć wszelkie wartości odstające w scenariuszu na żywo (za każdym razem, gdy rejestrowana jest nowa próbka). Jaki jest najlepszy sposób na osiągnięcie tych wyników?

Obecnie używam średniej ruchomej, aby usunąć trochę hałasu, ale co dalej? Proste rzeczy, takie jak odchylenie standardowe, szalenie, ... w stosunku do całego zestawu danych nie działa dobrze (nie mogę założyć, że szeregi czasowe są nieruchome) i chciałbym czegoś bardziej „dokładnego”, najlepiej czarnej skrzynki, takiej jak:

double outlier_detection (double * vector, double value);

gdzie wektor jest tablicą podwójnej zawierającej dane historyczne, a zwracana wartość jest wynikiem anomalii dla nowej „wartości” próbki.

Oto prosta funkcja R, która znajdzie wartości odstające szeregu czasowego (i opcjonalnie pokaże je na wykresie). Będzie obsługiwał sezonowe i nie sezonowe szeregi czasowe. Podstawową ideą jest znalezienie wiarygodnych oszacowań trendu i składników sezonowych i odjęcie ich. Następnie znajdź wartości odstające w resztkach. Test resztkowych wartości odstających jest taki sam, jak w przypadku standardowego wykresu pudełkowego - przyjmuje się, że punkty większe niż 1,5 IQR powyżej lub poniżej górnego i dolnego kwartylu są wartościami odstającymi. Liczba IQR powyżej / poniżej tych progów jest zwracana jako „wynik” odstający. Tak więc wynik może być dowolną liczbą dodatnią i będzie wynosił zero dla wartości nietypowych.

Zdaję sobie sprawę, że nie wdrażasz tego w języku R, ale często uważam, że funkcja R jest dobrym miejscem do rozpoczęcia. Następnie zadaniem jest przetłumaczenie tego na dowolny wymagany język.

tsoutliers <- function(x,plot=FALSE)
{
    x <- as.ts(x)
    if(frequency(x)>1)
        resid <- stl(x,s.window="periodic",robust=TRUE)$time.series[,3]
    else
    {
        tt <- 1:length(x)
        resid <- residuals(loess(x ~ tt))
    }
    resid.q <- quantile(resid,prob=c(0.25,0.75))
    iqr <- diff(resid.q)
    limits <- resid.q + 1.5*iqr*c(-1,1)
    score <- abs(pmin((resid-limits[1])/iqr,0) + pmax((resid - limits[2])/iqr,0))
    if(plot)
    {
        plot(x)
        x2 <- ts(rep(NA,length(x)))
        x2[score>0] <- x[score>0]
        tsp(x2) <- tsp(x)
        points(x2,pch=19,col="red")
        return(invisible(score))
    }
    else
        return(score)
}

Dobre rozwiązanie będzie miało kilka składników, w tym:

• Użyj odpornego, ruchomego okna gładkiego, aby usunąć niestabilność.

• Ponownie wyraż oryginalne dane, aby reszty względem gładkiego były w przybliżeniu symetrycznie rozmieszczone. Biorąc pod uwagę naturę twoich danych, prawdopodobne jest, że ich pierwiastki kwadratowe lub logarytmy dawałyby symetryczne reszty.

• Zastosuj metody kart kontrolnych lub przynajmniej myślenie kart kontrolnych do reszt.

O ile to ostatnie, myślenie na podstawie tabeli kontrolnej pokazuje, że „konwencjonalne” progi, takie jak 2 SD lub 1,5-krotność IQR poza kwartylami, działają słabo, ponieważ wyzwalają zbyt wiele fałszywych sygnałów wymykających się spod kontroli. Ludzie zwykle używają 3 SD w pracy na karcie kontrolnej, skąd 2,5 (lub nawet 3) razy IQR poza kwartylami byłoby dobrym punktem wyjścia.

Mniej więcej nakreśliłem naturę rozwiązania Roba Hyndmana, dodając do niego dwa główne punkty: potencjalną potrzebę ponownego wyrażenia danych i mądrość bycia bardziej konserwatywnym w sygnalizowaniu wartości odstającej. Nie jestem jednak pewien, czy Loess nadaje się do wykrywania online, ponieważ nie działa dobrze w punktach końcowych. Zamiast tego możesz użyć czegoś tak prostego jak ruchomy filtr środkowy (jak w przypadku wygładzania opornego Tukeya). Jeśli wartości odstające nie występują w seriach, możesz użyć wąskiego okna (być może 5 punktów danych, które rozpadnie się tylko w przypadku serii 3 lub więcej wartości odstających w grupie 5).

Po przeprowadzeniu analizy w celu ustalenia dobrego ponownego wyrażenia danych, jest mało prawdopodobne, że będziesz musiał zmienić to wyrażenie. Dlatego twój wykrywacz online naprawdę musi tylko odwoływać się do najnowszych wartości (ostatnie okno), ponieważ w ogóle nie użyje wcześniejszych danych. Jeśli masz naprawdę długie szeregi czasowe, możesz przejść dalej, aby przeanalizować autokorelację i sezonowość (takie jak powtarzające się codzienne lub tygodniowe wahania), aby usprawnić procedurę.

(Ta odpowiedź stanowi odpowiedź na zduplikowane (teraz zamknięte) pytanie w Wykrywanie zaległych zdarzeń , które przedstawiało niektóre dane w formie graficznej.)

Wykrywanie wartości odstających zależy od charakteru danych i tego, co chcesz o nich założyć. Metody ogólnego zastosowania opierają się na solidnych statystykach. Ideą tego podejścia jest scharakteryzowanie dużej ilości danych w sposób, na który nie mają wpływu żadne wartości odstające, a następnie wskazanie dowolnych indywidualnych wartości, które nie mieszczą się w tej charakterystyce.

Ponieważ jest to szereg czasowy, komplikuje to konieczność (ponownego) wykrywania wartości odstających na bieżąco. Jeśli ma to być zrobione, gdy seria się rozwija, to wolno nam wykorzystywać tylko starsze dane do wykrywania, a nie przyszłe dane! Ponadto, jako ochronę przed wieloma powtarzanymi testami, chcielibyśmy zastosować metodę, która ma bardzo niski odsetek wyników fałszywie dodatnich.

Te rozważania sugerują przeprowadzenie prostego, niezawodnego testu danych odstających od ruchomego okna na danych . Istnieje wiele możliwości, ale jedna prosta, łatwa do zrozumienia i łatwa do wdrożenia oparta jest na działającym MAD: absolutna mediana odchylenia od mediany. Jest to bardzo solidna miara zmienności danych, podobna do odchylenia standardowego. Źródło nadmiernych szczyt będzie kilka Mads lub bardziej większy niż mediana.

R$x = (1,2,\ldots,n)$$n=1150$$y$

# Parameters to tune to the circumstances:
window <- 30
threshold <- 5

# An upper threshold ("ut") calculation based on the MAD:
library(zoo) # rollapply()
ut <- function(x) {m = median(x); median(x) + threshold * median(abs(x - m))}
z <- rollapply(zoo(y), window, ut, align="right")
z <- c(rep(z[1], window-1), z) # Use z[1] throughout the initial period
outliers <- y > z

# Graph the data, show the ut() cutoffs, and mark the outliers:
plot(x, y, type="l", lwd=2, col="#E00000", ylim=c(0, 20000))
lines(x, z, col="Gray")
points(x[outliers], y[outliers], pch=19)

Zastosowany do zestawu danych, takiego jak czerwona krzywa zilustrowana w pytaniu, daje następujący wynik:

Dane są pokazane na czerwono, 30-dniowe okno mediany + 5 * progów MAD na szaro, a wartości odstające - które są po prostu wartościami danych powyżej szarej krzywej - na czarno.

(Próg można obliczyć dopiero na końcu początkowego okna. Dla wszystkich danych w tym początkowym oknie stosuje się pierwszy próg: dlatego szara krzywa jest płaska między x = 0 a x = 30).

Skutkami zmiany parametrów są (a) zwiększenie wartości windowbędzie miało tendencję do wygładzania szarej krzywej i (b) zwiększenie thresholdspowoduje podniesienie szarej krzywej. Wiedząc o tym, można wziąć początkowy segment danych i szybko zidentyfikować wartości parametrów, które najlepiej oddzielają odległe szczyty od reszty danych. Zastosuj te wartości parametrów, aby sprawdzić resztę danych. Jeśli wykres pokazuje, że metoda pogarsza się z czasem, oznacza to, że charakter danych się zmienia i parametry mogą wymagać ponownego dostrojenia.

Zauważ, jak mało ta metoda zakłada o danych: nie muszą one być normalnie dystrybuowane; nie muszą wykazywać żadnej częstotliwości; nie muszą nawet być nieujemne. Wszystko to zakłada się, że dane zachowują się w podobny sposób uzasadniony w czasie i że Dalekie szczyty są wyraźnie wyższe niż w pozostałej części danych.

Jeśli ktoś chciałby eksperymentować (lub porównać jakieś inne rozwiązanie z oferowanym tutaj), oto kod, którego użyłem do wygenerowania danych takich jak te pokazane w pytaniu.

n.length <- 1150
cycle.a <- 11
cycle.b <- 365/12
amp.a <- 800
amp.b <- 8000

set.seed(17)
x <- 1:n.length
baseline <- (1/2) * amp.a * (1 + sin(x * 2*pi / cycle.a)) * rgamma(n.length, 40, scale=1/40)
peaks <- rbinom(n.length, 1,  exp(2*(-1 + sin(((1 + x/2)^(1/5) / (1 + n.length/2)^(1/5))*x * 2*pi / cycle.b))*cycle.b))
y <- peaks * rgamma(n.length, 20, scale=amp.b/20) + baseline

Jeśli martwisz się założeniami z jakimś konkretnym podejściem, jednym z podejść jest wyszkolenie wielu uczniów na różnych sygnałach, a następnie użycie metod grupowania i zsumowanie „głosów” od twoich uczniów, aby dokonać klasyfikacji odstającej.

BTW, może to być warte przeczytania lub przeglądania, ponieważ odnosi się do kilku podejść do problemu.

• Wykrywanie wartości odstających online przez strumienie danych

Zgaduję, że wyrafinowany model szeregów czasowych nie będzie dla ciebie działał z powodu czasu potrzebnego na wykrycie wartości odstających przy użyciu tej metodologii. Dlatego oto obejście:

1. Najpierw określ bazowe „normalne” wzorce ruchu na rok na podstawie ręcznej analizy danych historycznych uwzględniających porę dnia, dzień tygodnia w weekend i miesiąc, rok itp.

2. Użyj tej linii bazowej wraz z prostym mechanizmem (np. Średnią ruchomą sugerowaną przez Carlosa), aby wykryć wartości odstające.

Możesz także przejrzeć literaturę dotyczącą statystycznego sterowania procesem w celu uzyskania pomysłów.

Sezonowo dostosowuj dane, aby normalny dzień wyglądał bardziej płasko. Możesz pobrać dzisiejszą próbkę o 17:00 i odjąć lub podzielić średnią z poprzednich 30 dni o 17:00. Następnie spójrz na odchylenia standardowe N (mierzone przy użyciu wstępnie dostosowanych danych) na wartości odstające. Można to zrobić osobno dla tygodniowych i codziennych „pór roku”.

Alternatywą dla podejścia nakreślonego przez Roba Hyndmana byłoby zastosowanie Prognozy Holta-Wintersa . Pasma ufności pochodzące z Holta-Wintersa można wykorzystać do wykrywania wartości odstających. Oto artykuł opisujący, jak używać Holt-Winters do „Wykrywania nieprawidłowych zachowań w szeregach czasowych do monitorowania sieci”. Implementacja dla RRDTool można znaleźć tutaj .

Analiza spektralna wykrywa okresowość w stacjonarnych szeregach czasowych. Podejście w dziedzinie częstotliwości oparte na estymacji gęstości widmowej to podejście, które poleciłbym jako pierwszy krok.

Jeżeli dla niektórych okresów nieregularność oznacza znacznie wyższy pik niż typowy dla tego okresu, wówczas szereg takich nieregularności nie byłby stacjonarny, a analiza spektralna nie byłaby odpowiednia. Zakładając, że określiłeś okres, w którym występują nieregularności, powinieneś być w stanie określić w przybliżeniu normalną wysokość piku, a następnie ustawić próg na pewnym poziomie powyżej tej średniej, aby wyznaczyć nieregularne przypadki.

Ponieważ są to dane szeregów czasowych, prosty filtr wykładniczy http://en.wikipedia.org/wiki/Exponential_smoothing wygładzi dane. Jest to bardzo dobry filtr, ponieważ nie trzeba gromadzić starych punktów danych. Porównaj każdą nowo wygładzoną wartość danych z jej niewygładzoną wartością. Gdy odchylenie przekroczy określony z góry próg (w zależności od tego, co uważasz za wartość odstającą w twoich danych), wtedy wartość odstająca może być łatwo wykryta.

W CI zrobi to dla próbki 16-bitowej w czasie rzeczywistym (uważam, że można ją znaleźć gdzieś tutaj <Objaśnienie - https://dsp.stackexchange.com/questions/378/what-is-the-best-first-order -iir-przybliżenie-do-ruchomej średniej filtru >)

#define BITS2 2     //< This is roughly = log2( 1 / alpha ), depending on how smooth you want your data to be

short Simple_Exp_Filter(int new_sample) 
{static int filtered_sample = 0;
long local_sample = sample << 16; /*We assume it is a 16 bit sample */
filtered_sample += (local_sample - filtered_sample) >> BITS2;   
return (short) ((filtered_sample+0x8000) >> 16); //< Round by adding .5 and truncating.   
}


int main()
{
newly_arrived = function_receive_new_sample();
filtered_sample = Simple_Exp_Filter(newly_arrived);
if (abs(newly_arrived - filtered_sample)/newly_arrived > THRESHOLD)
    {
    //AN OUTLIER HAS BEEN FOUND
    }
 return 0;   
}

Możesz użyć odchylenia standardowego ostatnich N pomiarów (musisz wybrać odpowiedni N). Dobrym wynikiem anomalii byłoby to, ile standardowych odchyleń stanowi pomiar od średniej ruchomej.

To, co robię, to grupowanie pomiarów według godziny i dnia tygodnia i porównanie standardowych odchyleń tego. Nadal nie poprawia takich rzeczy jak wakacje i sezonowość lato / zima, ale jest poprawna przez większość czasu.

Minusem jest to, że naprawdę musisz zebrać około roku danych, aby mieć wystarczająco dużo, aby stddev zaczął mieć sens.

Sugeruję poniższy schemat, który powinien być możliwy do wdrożenia za około jeden dzień:

Trening

• Zbierz jak najwięcej próbek, które możesz zachować w pamięci
• Usuń oczywiste wartości odstające, używając standardowego odchylenia dla każdego atrybutu
• Oblicz i przechowuj macierz korelacji, a także średnią każdego atrybutu
• Oblicz i przechowuj odległości Mahalanobisa wszystkich próbek

Obliczanie „odstających wartości”:

W przypadku pojedynczej próbki, dla której chcesz poznać jej „nietypowość”:

• Odzyskaj średnie, macierz kowariancji i dystans Mahalanobisa ze szkolenia
• Oblicz odległość Mahalanobisa „d” dla swojej próbki
• Zwraca percentyl, w którym spada „d” (używając odległości Mahalanobisa od treningu)

To będzie twój wynik odstający: 100% jest skrajnym wynikiem odstającym.

W przypadku, gdy trzeba szybko obliczyć wartości odstające, można skorzystać z pomysłu Roba Hyndmana i Mahito Sugiyamy ( https://github.com/BorgwardtLab/sampling-outlier-detection , biblioteka (spoutlier), funkcja qsp) do obliczenia wartości odstające, takie jak:

library(spoutlier)
rapidtsoutliers <- function(x,plot=FALSE,seed=123)
{
    set.seed(seed)
    x <- as.numeric(x)
    tt <- 1:length(x)
    qspscore <- qsp(x)
    limit <- quantile(qspscore,prob=c(0.95))
    score <- pmax((qspscore - limit),0)
    if(plot)
    {
        plot(x,type="l")
        x2 <- ts(rep(NA,length(x)))
        x2[score>0] <- x[score>0]
        tsp(x2) <- tsp(x)
        points(x2,pch=19,col="red")
        return(invisible(score))
    }
    else
        return(score)
}

wykrywanie anomalii wymaga budowy równania opisującego oczekiwanie. Wykrywanie interwencji jest dostępne zarówno w przypadku przyczynowym, jak i przyczynowym. Jeśli ktoś ma szereg prognostyczny, taki jak cena, sprawy mogą się nieco skomplikować. Inne odpowiedzi tutaj nie wydają się uwzględniać możliwej do przypisania przyczyny przypisywanej określonej przez użytkownika serii predyktorów, takiej jak cena, a zatem mogą być wadliwe. Ilość sprzedana może zależeć od ceny, być może wcześniejszych cen i być może ilości sprzedanej w przeszłości. Podstawę wykrycia anomalii (impulsy, impulsy sezonowe, zmiany poziomów i trendy czasu lokalnego) można znaleźć na stronie https://pdfs.semanticscholar.org/09c4/ba8dd3cc88289caf18d71e8985bdd11ad21c.pdf