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

16

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.

gianluca
quelle

Antworten:

6

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 (MM=1000Bi,ji=1Mj=1NNkBi,jj 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 .)(i1)/Mi/MN[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 Kjund dasjte Element der mittleren absoluten Abweichung wird geschätzt durch1

X¯j=1Kii1/2MBi,j,
j
1Ki|Xj¯i1/2M|Bi,j

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.

shabbychef
quelle
Wow, wirklich interessanter Ansatz. Das wusste ich nicht und ich werde es mir merken. Leider wird es in diesem Fall nicht funktionieren, da ich vom Standpunkt der Speichernutzung sehr restriktive Anforderungen habe, daher sollte M wirklich klein sein und ich denke, es würde zu viel Präzisionsverlust geben.
Gianluca
@gianluca: Es hört sich so an, als hätten Sie 1. viele Daten, 2. begrenzte Speicherressourcen, 3. hohe Genauigkeitsanforderungen. Ich kann sehen, warum dieses Problem Sie ausflippt! Vielleicht können Sie, wie von @kwak erwähnt, ein anderes Maß für die Streuung berechnen: MAD, IQR, Standardabweichung. Alle haben Ansätze, die für Ihr Problem funktionieren könnten.
Shabbychef
gianluca:> Geben Sie uns eine genauere Vorstellung von der Größe des Speichers, der Arrays und der Genauigkeit, die Sie wünschen. Es kann jedoch gut sein, dass Ihre Frage am besten über stackoverflow beantwortet wird.
User603
4

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.

fmark
quelle
2
Mir fehlt etwas: Wenn Sie "den gesamten Vektor im Speicher behalten" müssen, wie qualifiziert sich dies als "Online" -Algorithmus?
whuber
@whuber Nein, etwas fehlt nicht, ich denke, es ist kein Online-Algorithmus. Es benötigt O(n)Speicher und im schlimmsten Fall O (n) Zeit für jedes hinzugefügte Element. Bei normalverteilten Daten (und wahrscheinlich auch bei anderen Distributionen) funktioniert dies jedoch recht effizient.
Markieren Sie den
3

XXXss2/π

shabbychef
quelle
Das ist eine interessante Idee. Sie könnten es möglicherweise durch die Online-Erkennung von Ausreißern ergänzen und diese verwenden, um die Schätzung im Laufe der Zeit zu ändern.
whuber
Sie könnten wahrscheinlich die Methode von Welford verwenden, um die Standardabweichung online zu berechnen, die ich in meiner zweiten Antwort dokumentiert habe.
Markieren Sie den
1
Man sollte jedoch beachten, dass man auf diese Weise die Robustheit von Schätzern wie der expliziten MAD verlieren kann, die manchmal zu ihrer Wahl gegen einfachere Alternativen treiben.
Quarz,
2

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.

user603
quelle
Könnten Sie bitte die Beziehung zwischen MAD und den beiden Medianen detaillieren oder referenzieren?
Quarz,
Es ist wirklich die Formel des MAD: medich=1n|xich-medich=1n|(daher zwei Mediane)
user603
Hehm, ich meinte eigentlich, wenn Sie erklären können, wie diese Beziehung es den beiden Medianen erlaubt, gleichzeitig zu sein; diese scheinen mir abhängig zu sein, da sich die Eingaben zum äußeren Median bei jeder hinzugefügten Stichprobe zur inneren Berechnung ändern können. Wie würden Sie sie parallel durchführen?
Quarz
Ich muss auf die Binmedian-Arbeit zurückgreifen, um Einzelheiten zu erfahren ... aber einen berechneten Wert des Medians angeben (medich=1nxich) und einen neuen Wert von xn+1 Der Algorithmus könnte berechnen medich=1n+1xich viel schneller als Ö(n) durch Identifizieren des Fachs, in das xn+1gehört. Ich verstehe nicht, wie diese Einsicht in der verrückten Berechnung auf den äußeren Median verallgemeinert werden konnte.
user603
1

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.

historgram of differences

fmark
quelle
Dies gibt aus folgendem Grund keine genaue Antwort: ichxichichxich. Sie berechnen das erstere, während das OP das letztere wünscht.
Shabbychef
Ich stimme zu, dass die Methode ungenau ist. Ich bin jedoch mit Ihrer Diagnose der Ungenauigkeit nicht einverstanden. Welfords Methode zur Berechnung der Varianz, die nicht einmal ein Quadrat enthält, weist einen ähnlichen Fehler auf. Doch wie ngroß wird, der error/nbekommt verschwindend klein, überraschend schnell.
Markieren Sie den
Welfords Methode hat kein Quadrat, weil sie die Varianz berechnet, nicht die Standardabweichung. Es scheint, als würden Sie mit dem Quadrat die Standardabweichung schätzen, nicht die mittlere absolute Abweichung. Vermisse ich etwas?
Shabbychef
@shabbychef Bei jeder Welford-Iteration wird der Beitrag des neuen Datenpunkts zur absoluten Abweichung im Quadrat berechnet. Also nehme ich die Quadratwurzel jedes quadratischen Beitrags, um zur absoluten Abweichung zurückzukehren. Sie können beispielsweise feststellen, dass ich die Quadratwurzel des Deltas ziehe, bevor ich sie zur Abweichungssumme addiere, und nicht wie bei der Standardabweichung danach.
Markieren Sie den
3
Ich sehe das Problem; Welfords verschleiert das Problem mit dieser Methode: Die Online-Schätzung des Mittelwerts wird anstelle der endgültigen Schätzung des Mittelwerts verwendet. Während die Welford-Methode für die Varianz genau ist (bis zur Abrundung), ist dies bei dieser Methode nicht der Fall. Das Problem liegt nicht an der sqrtUngenauigkeit. Dies liegt daran, dass die Schätzung des laufenden Mittels verwendet wird. Um zu sehen, wann dies fehlschlägt, versuchen Sie Folgendes: xs <- sort(rnorm(n.testitems)) Wenn ich dies mit Ihrem Code versuche (nachdem ich ihn repariert habe, um zurückzukehren a.dev / n), erhalte ich relative Fehler in der Größenordnung von 9% -16%. Diese Methode ist also nicht permutationsinvariant, was zu Chaos führen könnte ...
shabbychef