Algorytm online dla średniego bezwzględnego odchylenia i dużego zestawu danych

Mam mały problem, który doprowadza mnie do szału. Muszę napisać procedurę dla procesu akwizycji online wielowymiarowych szeregów czasowych. Za każdym razem (na przykład 1 sekundę) otrzymuję nową próbkę, która jest w zasadzie wektorem zmiennoprzecinkowym o rozmiarze N. Operacja, którą muszę wykonać, jest nieco trudna:

1. Dla każdej nowej próbki obliczam wartości procentowe dla tej próbki (normalizując wektor, aby elementy sumowały się do 1).

2. Średni wektor procentowy obliczam w ten sam sposób, ale używając wcześniejszych wartości.

3. Dla każdej przeszłej wartości obliczam bezwzględne odchylenie wektora procentowego związane z tą próbką za pomocą globalnego średniego wektora procentowego obliczonego w kroku 2. W ten sposób absolutne odchylenie jest zawsze liczbą z zakresu od 0 (gdy wektor jest równy średniej wektor) i 2 (gdy jest całkowicie inny).

4. Używając średniej odchyleń dla wszystkich poprzednich próbek, obliczam średnie odchylenie bezwzględne, które ponownie jest liczbą między 0 a 2.

5. Używam średniego odchylenia bezwzględnego, aby wykryć, czy nowa próbka jest kompatybilna z innymi próbkami (porównując jej bezwzględne odchylenie ze średnim bezwzględnym odchyleniem całego zestawu obliczonego w kroku 4).

Ponieważ za każdym razem, gdy pobierana jest nowa próbka, zmiany globalnej średniej (a więc również zmiany średniej bezwzględnej odchyłki), czy istnieje sposób na obliczenie tej wartości bez wielokrotnego skanowania całego zestawu danych? (jeden raz do obliczenia średnich globalnych wartości procentowych i jeden raz do zebrania bezwzględnych odchyleń). Ok, wiem, że absolutnie łatwo jest obliczyć średnie globalne bez skanowania całego zestawu, ponieważ po prostu muszę użyć wektora tymczasowego do przechowywania sumy każdego wymiaru, ale co ze średnim absolutnym odchyleniem? Jego obliczenia obejmują abs()operatora, więc potrzebuję dostępu do wszystkich danych z przeszłości!

Dzięki za pomoc.

Jeśli możesz zaakceptować pewną niedokładność, problem ten można łatwo rozwiązać, licząc binning . Oznacza to, że piłki pewną liczbę largeish (powiedzmy, M = 1000 ), a następnie kilka pojemników zainicjować całkowite B I , j o i = 1 ... M , a j = 1 ... N , gdzie N jest wielkością wektora jako zero. Wtedy, gdy widzisz k th obserwację procentowy przyrost wektora, B í , j jeśli j th element tego wektora jest między $M$$M = 1000$$B_{i,j}$$i = 1\ldots M$$j = 1\ldots N$$N$$k$$B_{i,j}$$j$ i i /$(i-1)/M$$i/M$, zapętlając na $N$ elementów wektora. (Zakładam, że wektory wejściowe są nieujemne, więc kiedy obliczasz swoje „wartości procentowe”, wektory są w zakresie ).$[0,1]$

W dowolnym momencie możesz oszacować średni wektor z pojemników i średnie bezwzględne odchylenie. Po stwierdzeniu, takich wektorów The J ty element średniej szacuje ˉ X J = $K$$j$iJelementem Th średnią bezwzględną odchylenia szacuje

edycja : jest to szczególny przypadek bardziej ogólnego podejścia, w którym buduje się szacunkową gęstość empiryczną. Można to zrobić za pomocą wielomianów, splajnów itp., Ale podejście grupowania jest najłatwiejsze do opisania i wdrożenia.

W przeszłości stosowałem następujące podejście do umiarkowanie wydajnego obliczania odchylenia rozgrzeszenia (zauważ, że jest to podejście programistów, a nie statystyków, więc niewątpliwie mogą istnieć sprytne sztuczki, takie jak shabbychef, które mogą być bardziej wydajne).

OSTRZEŻENIE: To nie jest algorytm online. Wymaga O(n)pamięci. Co więcej, ma najgorszy wynik w O(n)przypadku takich zestawów danych [1, -2, 4, -8, 16, -32, ...](tj. Taki sam jak w przypadku pełnego przeliczenia). [1]

Ponieważ jednak nadal działa dobrze w wielu przypadkach użycia, warto opublikować tutaj. Na przykład, aby obliczyć absolutne odchylenie 10000 losowo liczb od -100 do 100 w miarę dostarczania każdego elementu, mój algorytm zajmuje mniej niż jedną sekundę, podczas gdy pełne ponowne obliczenie zajmuje ponad 17 sekund (na mojej maszynie będą się różnić w zależności od maszyny i zgodnie z danymi wejściowymi). Musisz jednak zachować cały wektor w pamięci, co może być ograniczeniem dla niektórych zastosowań. Zarys algorytmu jest następujący:

1. Zamiast pojedynczego wektora do przechowywania poprzednich pomiarów, użyj trzech posortowanych kolejek priorytetowych (coś w rodzaju sterty min / max). Te trzy listy dzielą dane wejściowe na trzy: pozycje większe niż średnia, pozycje mniejsze niż średnia i pozycje równe średniej.
2. (Prawie) za każdym razem, gdy dodajesz przedmiot, zmienia się średnia, więc musimy podzielić na części. Kluczową sprawą jest posortowana natura partycji, co oznacza, że ​​zamiast skanować każdy element na liście do podziału, musimy tylko czytać te elementy, które przenosimy. W najgorszym przypadku będzie to nadal wymagaćO(n) operacji przenoszenia, w wielu przypadkach użycia tak nie jest.
3. Używając sprytnej księgowości, możemy upewnić się, że odchylenie jest poprawnie obliczane przez cały czas, podczas podziału partycji i dodawania nowych elementów.

Przykładowy kod w pythonie znajduje się poniżej. Pamiętaj, że pozwala tylko dodawać elementy do listy, a nie usuwać. Można to łatwo dodać, ale w chwili, gdy to pisałem, nie było takiej potrzeby. Zamiast samodzielnie wdrażać kolejki priorytetowe, skorzystałem z sortowanej listy z doskonałego pakietu blist Daniela Stutzbacha , który wykorzystuje B + Tree .

Rozważ ten kod na licencji MIT . Nie został znacznie zoptymalizowany ani dopracowany, ale działał dla mnie w przeszłości. Nowe wersje będą dostępne tutaj . Daj mi znać, jeśli masz jakieś pytania lub znajdziesz jakieś błędy.

from blist import sortedlist
import operator

class deviance_list:
    def __init__(self):
        self.mean =  0.0
        self._old_mean = 0.0
        self._sum =  0L
        self._n =  0  #n items
        # items greater than the mean
        self._toplist =  sortedlist()
        # items less than the mean
        self._bottomlist = sortedlist(key = operator.neg)
        # Since all items in the "eq list" have the same value (self.mean) we don't need
        # to maintain an eq list, only a count
        self._eqlistlen = 0

        self._top_deviance =  0
        self._bottom_deviance =  0

    @property
    def absolute_deviance(self):
        return self._top_deviance + self._bottom_deviance

    def append(self,  n):
        # Update summary stats
        self._sum += n
        self._n +=  1
        self._old_mean =  self.mean
        self.mean =  self._sum /  float(self._n)

        # Move existing things around
        going_up = self.mean > self._old_mean
        self._rebalance(going_up)

        # Add new item to appropriate list
        if n >  self.mean:
            self._toplist.add(n)
            self._top_deviance +=  n -  self.mean
        elif n == self.mean: 
            self._eqlistlen += 1
        else:
            self._bottomlist.add(n)
            self._bottom_deviance += self.mean -  n


    def _move_eqs(self,  going_up):
        if going_up:
            self._bottomlist.update([self._old_mean] *  self._eqlistlen)
            self._bottom_deviance += (self.mean - self._old_mean) * self._eqlistlen
            self._eqlistlen = 0
        else:
            self._toplist.update([self._old_mean] *  self._eqlistlen)
            self._top_deviance += (self._old_mean - self.mean) * self._eqlistlen
            self._eqlistlen = 0


    def _rebalance(self, going_up):
        move_count,  eq_move_count = 0, 0
        if going_up:
            # increase the bottom deviance of the items already in the bottomlist
            if self.mean !=  self._old_mean:
                self._bottom_deviance += len(self._bottomlist) *  (self.mean -  self._old_mean)
                self._move_eqs(going_up)


            # transfer items from top to bottom (or eq) list, and change the deviances
            for n in iter(self._toplist):
                if n < self.mean:
                    self._top_deviance -= n -  self._old_mean
                    self._bottom_deviance += (self.mean -  n)
                    # we increment movecount and move them after the list
                    # has finished iterating so we don't modify the list during iteration
                    move_count +=  1
                elif n == self.mean:
                    self._top_deviance -= n -  self._old_mean
                    self._eqlistlen += 1
                    eq_move_count +=  1
                else:
                    break
            for _ in xrange(0,  move_count):
                self._bottomlist.add(self._toplist.pop(0))
            for _ in xrange(0,  eq_move_count):
                self._toplist.pop(0)

            # decrease the top deviance of the items remain in the toplist
            self._top_deviance -= len(self._toplist) *  (self.mean -  self._old_mean)
        else:
            if self.mean !=  self._old_mean:
                self._top_deviance += len(self._toplist) *  (self._old_mean -  self.mean)
                self._move_eqs(going_up)
            for n in iter(self._bottomlist): 
                if n > self.mean:
                    self._bottom_deviance -= self._old_mean -  n
                    self._top_deviance += n -  self.mean
                    move_count += 1
                elif n == self.mean:
                    self._bottom_deviance -= self._old_mean -  n
                    self._eqlistlen += 1
                    eq_move_count +=  1
                else:
                    break
            for _ in xrange(0,  move_count):
                    self._toplist.add(self._bottomlist.pop(0))
            for _ in xrange(0,  eq_move_count):
                self._bottomlist.pop(0)

            # decrease the bottom deviance of the items remain in the bottomlist
            self._bottom_deviance -= len(self._bottomlist) *  (self._old_mean -  self.mean)


if __name__ ==  "__main__":
    import random
    dv =  deviance_list()
    # Test against some random data,  and calculate result manually (nb. slowly) to ensure correctness
    rands = [random.randint(-100,  100) for _ in range(0,  1000)]
    ns = []
    for n in rands: 
        dv.append(n)
        ns.append(n)
        print("added:%4d,  mean:%3.2f,  oldmean:%3.2f,  mean ad:%3.2f" %
              (n, dv.mean,  dv._old_mean,  dv.absolute_deviance / dv.mean))
        assert sum(ns) == dv._sum,  "Sums not equal!"
        assert len(ns) == dv._n,  "Counts not equal!"
        m = sum(ns) / float(len(ns))
        assert m == dv.mean,  "Means not equal!"
        real_abs_dev = sum([abs(m - x) for x in ns])
        # Due to floating point imprecision, we check if the difference between the
        # two ways of calculating the asb. dev. is small rather than checking equality
        assert abs(real_abs_dev - dv.absolute_deviance) < 0.01, (
            "Absolute deviances not equal. Real:%.2f,  calc:%.2f" %  (real_abs_dev,  dv.absolute_deviance))

[1] Jeśli objawy utrzymują się, skontaktuj się z lekarzem.

Istnieje również podejście parametryczne. Ignorując wektorową naturę danych i patrząc tylko na marginesy, wystarczy rozwiązać problem: znajdź internetowy algorytm do obliczenia średniego bezwzględnego odchylenia skalara$X$. Jeśli (i to jest tutaj duże „jeśli”), tak myślałeś$X$po pewnym rozkładzie prawdopodobieństwa o nieznanych parametrach można oszacować parametry za pomocą algorytmu online, a następnie obliczyć średnie bezwzględne odchylenie na podstawie tego sparametryzowanego rozkładu. Na przykład, jeśli tak myślałeś$X$ był (w przybliżeniu) normalnie rozłożony, można oszacować jego odchylenie standardowe, jak $s$, a średnie bezwzględne odchylenie zostanie oszacowane na podstawie $s \sqrt{2 / \pi}$(patrz Połowa rozkładu normalnego ).

MAD (x) to tylko dwa równoległe obliczenia mediany, z których każde można wykonać online za pomocą algorytmu binmedian .

Powiązany artykuł, a także kod C i FORTRAN można znaleźć tutaj .

(jest to po prostu zastosowanie sprytnej sztuczki na szczycie sprytnej sztuczki Shabbychef, aby zaoszczędzić na pamięci).

Uzupełnienie:

Istnieje wiele starszych wieloprzebiegowych metod obliczania kwantyli. Popularnym podejściem jest utrzymywanie / aktualizowanie wyznaczonego rozmiaru zbiornika obserwacji losowo wybranych ze strumienia i rekurencyjne obliczanie kwantyli (patrz ten przegląd) na tym zbiorniku. To (i powiązane) podejście zastępuje to zaproponowane powyżej.

Poniżej podano niedokładne przybliżenie, chociaż niedokładność będzie zależeć od rozkładu danych wejściowych. Jest to algorytm online, ale tylko przybliża absolutne odchylenie. Opiera się na dobrze znanym algorytmie obliczania wariancji online, opisanym przez Welforda w latach 60. Jego algorytm, przetłumaczony na R, wygląda następująco:

M2 <- 0
mean <- 0
n <- 0

var.online <- function(x){
    n <<- n + 1
    diff <- x - mean
    mean <<- mean + diff / n
    M2 <<- M2 + diff * (x - mean)
    variance <- M2 / (n - 1)
    return(variance)
}

Działa bardzo podobnie do wbudowanej funkcji wariancji R:

set.seed(2099)
n.testitems <- 1000
n.tests <- 100
differences <- rep(NA, n.tests)
for (i in 1:n.tests){
        # Reset counters
        M2 <- 0
        mean <- 0
        n <- 0

        xs <- rnorm(n.testitems)
        for (j in 1:n.testitems){
                v <- var.online(xs[j])
        }

        differences[i] <- abs(v - var(xs))

}
summary(differences)
     Min.   1st Qu.    Median      Mean   3rd Qu.      Max. 
0.000e+00 2.220e-16 4.996e-16 6.595e-16 9.992e-16 1.887e-15 

Modyfikacja algorytmu w celu obliczenia absolutnego odchylenia wymaga po prostu dodatkowego sqrtwywołania. Jednakże sqrtwprowadza nieścisłości, które znajdują odzwierciedlenie w wyniku:

absolute.deviance.online <- function(x){
    n <<- n + 1
    diff <- x - mean
    mean <<- mean + diff / n
    a.dev <<- a.dev + sqrt(diff * (x - mean))
    return(a.dev)
}

Błędy, obliczone jak wyżej, są znacznie większe niż w przypadku obliczania wariancji:

    Min.  1st Qu.   Median     Mean  3rd Qu.     Max. 
0.005126 0.364600 0.808000 0.958800 1.360000 3.312000 

Jednak w zależności od przypadku użycia ta wielkość błędu może być do zaakceptowania.