Online-Algorithmus für mittlere absolute Abweichung und großen Datensatz

Ich habe ein kleines Problem, das mich ausflippen lässt. Ich muss eine Prozedur für einen Online-Erfassungsprozess einer multivariaten Zeitreihe schreiben. In jedem Zeitintervall (zum Beispiel 1 Sekunde) erhalte ich eine neue Stichprobe, die im Grunde genommen ein Gleitkomma-Vektor der Größe N ist. Die Operation, die ich ausführen muss, ist etwas schwierig:

1. Für jede neue Stichprobe berechne ich die Prozentsätze für diese Stichprobe (indem ich den Vektor so normalisiere, dass die Elemente zu 1 summieren).

2. Ich berechne den durchschnittlichen prozentualen Vektor auf die gleiche Weise, aber unter Verwendung der vergangenen Werte.

3. Für jeden vergangenen Wert berechne ich die absolute Abweichung des mit dieser Stichprobe verbundenen Prozentvektors mit dem in Schritt 2 berechneten globalen durchschnittlichen Prozentvektor. Auf diese Weise ist die absolute Abweichung immer eine Zahl zwischen 0 (wenn der Vektor dem Durchschnitt entspricht) Vektor) und 2 (wenn es total verschieden ist).

4. Aus dem Durchschnitt der Abweichungen aller vorhergehenden Stichproben berechne ich die mittlere absolute Abweichung, die wiederum eine Zahl zwischen 0 und 2 ist.

5. Ich verwende die mittlere absolute Abweichung, um festzustellen, ob eine neue Probe mit den anderen Proben kompatibel ist (indem ich ihre absolute Abweichung mit der mittleren absoluten Abweichung des gesamten in Schritt 4 berechneten Satzes vergleiche).

Gibt es eine Möglichkeit, diesen Wert zu berechnen, ohne den gesamten Datensatz mehrmals zu scannen, da sich jedes Mal, wenn eine neue Stichprobe erfasst wird, der globale Durchschnitt ändert (und sich somit auch die mittlere absolute Abweichung ändert)? (einmal für die Berechnung der globalen Durchschnittsprozentsätze und einmal für die Erfassung der absoluten Abweichungen). Ok, ich weiß, dass es absolut einfach ist, die globalen Mittelwerte zu berechnen, ohne den gesamten Satz abzutasten, da ich nur einen temporären Vektor zum Speichern der Summe jeder Dimension verwenden muss, aber was ist mit der mittleren absoluten Abweichung? abs()Da die Berechnung den Operator einschließt , muss ich auf alle früheren Daten zugreifen können!

Danke für Ihre Hilfe.

Wenn Sie eine gewisse Ungenauigkeit akzeptieren, kann dieses Problem leicht durch Binning- Zählungen gelöst werden . Das heißt, wählen Sie eine größere Zahl (sagen wir M = 1000 ) und initialisieren Sie dann einige ganzzahlige Bins B i , j für i = 1 … M und j = 1 … N , wobei N die Vektorgröße ist, als Null. Wenn Sie dann die k- te Beobachtung eines prozentualen Vektors sehen, erhöhen Sie B i , j, wenn das j- te Element dieses Vektors zwischen $M$$M = 1000$$B_{i,j}$$i = 1\ldots M$$j = 1\ldots N$$N$$k$$B_{i,j}$$j$ und i / M , wobei die N Elemente des Vektors durchlaufen werden. (Ich gehe davon aus, dass Ihre Eingabevektoren nicht negativ sind, sodass die Vektoren bei der Berechnung Ihrer 'Prozentzahlen' im Bereich [ 0 , 1 ] liegen .)$(i-1)/M$$i/M$$N$$[0,1]$

Zu jedem Zeitpunkt können Sie den mittleren Vektor aus den Behältern und die mittlere absolute Abweichung schätzen. Nach der Beobachtung von solchen Vektoren wird das j- te Element des Mittelwerts mit ˉ X j = 1 geschätzt $K$$j$und dasjte Element der mittleren absoluten Abweichung wird geschätzt durch

Bearbeiten : Dies ist ein spezieller Fall eines allgemeineren Ansatzes, bei dem Sie eine empirische Dichteschätzung erstellen. Dies könnte mit Polynomen, Splines usw. geschehen, aber der Binning-Ansatz ist am einfachsten zu beschreiben und zu implementieren.

Ich habe in der Vergangenheit den folgenden Ansatz verwendet, um die Absolutionsabweichung mäßig effizient zu berechnen (beachten Sie, dass dies ein Programmieransatz ist, kein Statistiker, so dass es zweifellos clevere Tricks wie Shabbychefs gibt , die effizienter sein könnten).

WARNUNG: Dies ist kein Online-Algorithmus. Es erfordert O(n)Speicher. Darüber hinaus hat es eine Worst-Case-Performance von O(n), für Datensätze wie [1, -2, 4, -8, 16, -32, ...](dh die gleiche wie die vollständige Neuberechnung). [1]

Da es jedoch in vielen Anwendungsfällen immer noch eine gute Leistung bringt, kann es sich lohnen, es hier zu veröffentlichen. Um zum Beispiel die absolute Abweichung von 10000 Zufallszahlen zwischen -100 und 100 zu berechnen , benötigt mein Algorithmus weniger als eine Sekunde, während die vollständige Neuberechnung über 17 Sekunden dauert (auf meinem Computer variieren je nach Computer und entsprechend den eingegebenen Daten). Sie müssen jedoch den gesamten Vektor im Speicher behalten, was für einige Anwendungen eine Einschränkung sein kann. Der Umriss des Algorithmus lautet wie folgt:

1. Verwenden Sie statt eines einzelnen Vektors zum Speichern vergangener Messungen drei sortierte Prioritätswarteschlangen (so etwas wie einen Min / Max-Heap). Diese drei Listen unterteilen die Eingabe in drei Elemente: Elemente, die größer als der Mittelwert sind, Elemente, die kleiner als der Mittelwert sind, und Elemente, die dem Mittelwert entsprechen.
2. (Fast) jedes Mal, wenn Sie ein Element hinzufügen, ändert sich der Mittelwert, sodass wir neu partitionieren müssen. Das Entscheidende ist die Sortierung der Partitionen. Dies bedeutet, dass wir nicht jedes Element in der Liste nach Partitionen durchsuchen, sondern nur die Elemente lesen müssen, die wir gerade verschieben. Während dies im schlimmsten Fall immer noch O(n)Verschiebevorgänge erfordert, ist dies in vielen Anwendungsfällen nicht der Fall.
3. Mit einer ausgeklügelten Buchführung können wir sicherstellen, dass die Abweichung jederzeit korrekt berechnet wird, wenn wir neu partitionieren und neue Elemente hinzufügen.

Einige Beispielcodes in Python finden Sie weiter unten. Beachten Sie, dass nur Elemente zur Liste hinzugefügt und nicht entfernt werden können. Dies könnte leicht hinzugefügt werden, aber zu dem Zeitpunkt, als ich dies schrieb, hatte ich keine Notwendigkeit dafür. Anstatt die Priority Queues selbst zu implementieren, habe ich die sortierte Liste aus Daniel Stutzbachs exzellentem Blist-Paket verwendet , die intern B + Tree s verwendet.

Betrachten Sie diesen Code als unter der MIT-Lizenz lizenziert . Es wurde nicht wesentlich optimiert oder poliert, hat aber in der Vergangenheit für mich funktioniert. Neue Versionen verfügbar sein hier . Lassen Sie mich wissen, wenn Sie Fragen haben oder Fehler finden.

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] Wenn die Symptome anhalten, wenden Sie sich an Ihren Arzt.

$X$$X$$X$$s$$s \sqrt{2 / \pi}$

MAD (x) ist nur zwei gleichzeitige Median-Berechnungen, von denen jede online über den Binmedian- Algorithmus durchgeführt werden kann.

Den dazugehörigen Artikel sowie den C- und FORTRAN-Code finden Sie hier online .

(Dies ist nur die Verwendung eines cleveren Tricks zusätzlich zu Shabbychefs cleverem Trick, um Speicherplatz zu sparen).

Nachtrag:

Es gibt eine Vielzahl älterer Multi-Pass-Methoden zur Berechnung von Quantilen. Ein gängiger Ansatz besteht darin, ein deterministisch großes Reservoir von Beobachtungen zu verwalten / zu aktualisieren, die zufällig aus dem Datenstrom ausgewählt wurden, und Quantile (siehe diese Übersicht) auf diesem Reservoir rekursiv zu berechnen . Dieser (und verwandte) Ansatz wird durch den oben vorgeschlagenen ersetzt.

Das Folgende liefert eine ungenaue Annäherung, obwohl die Ungenauigkeit von der Verteilung der Eingabedaten abhängt. Es ist ein Online-Algorithmus, der sich jedoch nur der absoluten Abweichung annähert. Es basiert auf einem bekannten Algorithmus zur Online- Varianzberechnung , der von Welford in den 1960er Jahren beschrieben wurde. Sein in R übersetzter Algorithmus sieht folgendermaßen aus:

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

Es funktioniert sehr ähnlich wie die in R integrierte Varianzfunktion:

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 

Das Ändern des Algorithmus zur Berechnung der absoluten Abweichung erfordert lediglich einen zusätzlichen sqrtAufruf. Diesqrt führt jedoch zu Ungenauigkeiten, die sich im Ergebnis widerspiegeln:

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

Die wie oben berechneten Fehler sind viel größer als bei der Varianzberechnung:

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

Abhängig von Ihrem Anwendungsfall kann diese Fehlergröße jedoch akzeptabel sein.