Algorithme en ligne pour l'écart absolu moyen et un grand ensemble de données

J'ai un petit problème qui me fait paniquer. Je dois écrire une procédure pour un processus d'acquisition en ligne d'une série temporelle multivariée. À chaque intervalle de temps (par exemple 1 seconde), j'obtiens un nouvel échantillon, qui est essentiellement un vecteur à virgule flottante de taille N. L'opération que je dois faire est un peu délicate:

1. Pour chaque nouvel échantillon, je calcule les pourcentages pour cet échantillon (en normalisant le vecteur pour que les éléments totalisent 1).

2. Je calcule le vecteur des pourcentages moyens de la même manière, mais en utilisant les valeurs passées.

3. Pour chaque valeur passée, je calcule l'écart absolu du vecteur de pourcentages lié à cet échantillon avec le vecteur de pourcentages moyen global calculé à l'étape 2. De cette façon, l'écart absolu est toujours un nombre compris entre 0 (lorsque le vecteur est égal à la moyenne vecteur) et 2 (quand il est totalement différent).

4. En utilisant la moyenne des écarts pour tous les échantillons précédents, je calcule l'écart absolu moyen, qui est à nouveau un nombre compris entre 0 et 2.

5. J'utilise l'écart absolu moyen pour détecter si un nouvel échantillon est compatible avec les autres échantillons (en comparant son écart absolu avec l'écart absolu moyen de l'ensemble calculé à l'étape 4).

Étant donné que chaque fois qu'un nouvel échantillon est collecté, les variations moyennes globales (et donc l'écart absolu moyen change également), existe-t-il un moyen de calculer cette valeur sans analyser plusieurs fois l'ensemble des données? (une fois pour le calcul des pourcentages moyens mondiaux et une fois pour la collecte des écarts absolus). Ok, je sais qu'il est absolument facile de calculer les moyennes globales sans balayer l'ensemble, car je n'ai qu'à utiliser un vecteur temporaire pour stocker la somme de chaque dimension, mais qu'en est-il de l'écart absolu moyen? Son calcul inclut l' abs()opérateur, j'ai donc besoin d'accéder à toutes les données passées!

Merci de votre aide.

Si vous pouvez accepter une certaine imprécision, ce problème peut être résolu facilement en binning compte. C'est-à-dire, choisissez un grand nombre (disons M = 1000 ), puis initialisez quelques bacs entiers B i , j pour i = 1 … M et j = 1 … N , où N est la taille du vecteur, comme zéro. Ensuite, lorsque vous voyez la k ème observation d'un vecteur en pourcentage, incrémentez B i , j si le j ème élément de ce vecteur est compris entre $M$$M = 1000$$B_{i,j}$$i = 1\ldots M$$j = 1\ldots N$$N$$k$$B_{i,j}$$j$ et i / M , en boucle sur les N éléments du vecteur. (Je suppose que vos vecteurs d'entrée ne sont pas négatifs, de sorte que lorsque vous calculez vos «pourcentages», les vecteurs sont dans la plage [ 0 , 1 ] .)$(i-1)/M$$i/M$$N$$[0,1]$

À tout moment, vous pouvez estimer le vecteur moyen des bacs et l'écart absolu moyen. Après avoir observé tels vecteurs, le j ème élément de la moyenne est estimé par ˉ X j = $K$$j$etjélément e de l'écart absolu moyen est estimé par

modifier : il s'agit d'un cas spécifique d'une approche plus générale où vous construisez une estimation de densité empirique. Cela pourrait être fait avec des polynômes, des splines, etc., mais l'approche de regroupement est la plus simple à décrire et à mettre en œuvre.

J'ai utilisé l'approche suivante dans le passé pour calculer la déviation d'absolution de manière modérément efficace (notez que c'est une approche de programmeurs, pas de statisticiens, donc il est indubitable qu'il peut y avoir des astuces intelligentes comme celles de shabbychef qui pourraient être plus efficaces).

AVERTISSEMENT: ce n'est pas un algorithme en ligne. Cela nécessite de la O(n)mémoire. De plus, il présente les performances les plus défavorables de O(n), pour des ensembles de données comme [1, -2, 4, -8, 16, -32, ...](c'est-à-dire les mêmes que pour le recalcul complet). [1]

Cependant, comme il fonctionne toujours bien dans de nombreux cas d'utilisation, il peut être utile de le publier ici. Par exemple, afin de calculer la déviance absolue de 10000 nombres aléatoires entre -100 et 100 à l'arrivée de chaque élément, mon algorithme prend moins d'une seconde, tandis que le recalcul complet prend plus de 17 secondes (sur ma machine, variera selon la machine et selon les données d'entrée). Cependant, vous devez conserver le vecteur entier en mémoire, ce qui peut être une contrainte pour certaines utilisations. Le contour de l'algorithme est le suivant:

1. Au lieu d'avoir un seul vecteur pour stocker les mesures passées, utilisez trois files d'attente prioritaires triées (quelque chose comme un tas min / max). Ces trois listes divisent l'entrée en trois: éléments supérieurs à la moyenne, éléments inférieurs à la moyenne et éléments égaux à la moyenne.
2. (Presque) chaque fois que vous ajoutez un élément, la moyenne change, nous devons donc repartitionner. La chose cruciale est la nature triée des partitions, ce qui signifie qu'au lieu d'analyser chaque élément de la liste pour les répartir, nous n'avons qu'à lire les éléments que nous déplaçons. Alors que dans le pire des cas, cela nécessitera toujours des O(n)opérations de déplacement, pour de nombreux cas d'utilisation, ce n'est pas le cas.
3. En utilisant une comptabilité intelligente, nous pouvons nous assurer que la déviance est correctement calculée à tout moment, lors du repartitionnement et lors de l'ajout de nouveaux éléments.

Un exemple de code, en python, est ci-dessous. Notez qu'il permet uniquement d'ajouter des éléments à la liste, pas de les supprimer. Cela pourrait facilement être ajouté, mais au moment où j'ai écrit cela, je n'en avais pas besoin. Plutôt que d'implémenter moi-même les files d'attente prioritaires, j'ai utilisé la liste de tri de l'excellent package blist de Daniel Stutzbach , qui utilise B + Tree s en interne.

Considérez ce code sous licence MIT . Il n'a pas été optimisé ou poli de manière significative, mais a fonctionné pour moi dans le passé. De nouvelles versions seront disponibles ici . Faites-moi savoir si vous avez des questions ou si vous trouvez des bugs.

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] Si les symptômes persistent, consultez votre médecin.

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

MAD (x) n'est que deux calculs médians simultanés, chacun pouvant être effectué en ligne via l' algorithme binmédian .

Vous pouvez trouver le papier associé ainsi que le code C et FORTRAN en ligne ici .

(c'est juste l'utilisation d'une astuce intelligente en plus de l'astuce intelligente de Shabbychef, pour économiser de la mémoire).

Addenda:

Il existe une multitude d'anciennes méthodes multi-passes pour le calcul des quantiles. Une approche populaire consiste à maintenir / mettre à jour un réservoir d'observations de taille déterministe sélectionné au hasard dans le cours d'eau et à calculer récursivement des quantiles (voir cette revue) sur ce réservoir. Cette approche (et connexe) est remplacée par celle proposée ci-dessus.

The following provides an inaccurate approximation, although the inaccuracy will depend on the distribution of the input data. It is an online algorithm, but only approximates the absolute deviance. It is based on a well known algorithm for calculating variance online, described by Welford in the 1960s. His algorithm, translated into R, looks like:

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

It performs very similarly to R's builtin variance function:

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 

Modifying the algorithm to calculate absolute deviation simply involves an additional sqrt call. However, the sqrt introduces inaccuracies that are reflected in the result:

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

The errors, calculated as above, are much greater than for the variance calculation:

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

However, depending on your use case, this magnitude of error might be acceptable.