Algoritmo en línea para desviación media absoluta y gran conjunto de datos

Tengo un pequeño problema que me está volviendo loco. Tengo que escribir el procedimiento para un proceso de adquisición en línea de una serie de tiempo multivariante. En cada intervalo de tiempo (por ejemplo, 1 segundo), obtengo una nueva muestra, que es básicamente un vector de coma flotante de tamaño N. La operación que necesito hacer es un poco complicada:

1. Para cada nueva muestra, calculo los porcentajes para esa muestra (normalizando el vector para que los elementos sumen 1).

2. Calculo el vector de porcentajes promedio de la misma manera, pero usando los valores pasados.

3. Para cada valor pasado, calculo la desviación absoluta del vector de porcentajes relacionado con esa muestra con el vector de porcentajes promedio global calculado en el paso 2. De esta manera, la desviación absoluta es siempre un número entre 0 (cuando el vector es igual al promedio vector) y 2 (cuando es totalmente diferente).

4. Utilizando el promedio de las desviaciones para todas las muestras anteriores, calculo la desviación absoluta media, que nuevamente es un número entre 0 y 2.

5. Utilizo la desviación media absoluta para detectar si una nueva muestra es compatible con las otras muestras (comparando su desviación absoluta con la desviación media absoluta de todo el conjunto calculado en el paso 4).

Dado que cada vez que se recolecta una nueva muestra, el promedio global cambia (y también la desviación absoluta media también cambia), ¿hay alguna manera de calcular este valor sin escanear el conjunto de datos completo varias veces? (una vez para calcular los porcentajes promedio globales y una vez para recopilar las desviaciones absolutas). Ok, sé que es absolutamente fácil calcular los promedios globales sin escanear todo el conjunto, ya que solo tengo que usar un vector temporal para almacenar la suma de cada dimensión, pero ¿qué pasa con la desviación absoluta media? Su cálculo incluye elabs() operador, ¡así que necesito acceder a todos los datos pasados!

Gracias por tu ayuda.

Si puede aceptar alguna imprecisión, este problema se puede resolver fácilmente mediante recuentos de agrupamiento . Es decir, elija un número grande (digamos M = 1000 ), luego inicialice algunos contenedores enteros B i , j para i = 1 ... M y j = 1 ... N , donde N es el tamaño del vector, como cero. Luego, cuando se ve la k ésima observación de un vector porcentual, incremento B i , j si el j -ésimo elemento de este vector es entre $M$$M = 1000$$B_{i,j}$$i = 1\ldots M$$j = 1\ldots N$$N$$k$$B_{i,j}$$j$). e i / M , recorriendo los N elementos del vector. (Supongo que sus vectores de entrada no son negativos, de modo que cuando calcula sus 'porcentajes', los vectores están en el rango [ 0 , 1 $(i-1)/M$$i/M$$N$$[0,1]$

En cualquier momento, puede estimar el vector medio de los contenedores y la desviación absoluta media. Después de observar tales vectores, el elemento j de la media se estima por ˉ X j = $K$$j$y elelementojde la desviación absoluta media se estima en

editar : este es un caso específico de un enfoque más general en el que está construyendo una estimación de densidad empírica. Esto podría hacerse con polinomios, splines, etc., pero el enfoque de agrupamiento es el más fácil de describir e implementar.

He utilizado el siguiente enfoque en el pasado para calcular la desviación de absolución de manera moderadamente eficiente (tenga en cuenta que este es un enfoque de programadores, no un estadístico, por lo que indudablemente puede haber trucos inteligentes como el de shabbychef que podrían ser más eficientes).

ADVERTENCIA: este no es un algoritmo en línea. Requiere O(n)memoria Además, tiene un rendimiento en el peor de los casos O(n), para conjuntos de datos como [1, -2, 4, -8, 16, -32, ...](es decir, el mismo que el recálculo completo). [1]

Sin embargo, debido a que aún funciona bien en muchos casos de uso, podría valer la pena publicarlo aquí. Por ejemplo, para calcular la desviación absoluta de 10000 números aleatorios entre -100 y 100 a medida que llega cada elemento, mi algoritmo tarda menos de un segundo, mientras que el recálculo completo lleva más de 17 segundos (en mi máquina, variará por máquina y según datos de entrada). Sin embargo, debe mantener todo el vector en la memoria, lo que puede ser una restricción para algunos usos. El esquema del algoritmo es el siguiente:

1. En lugar de tener un solo vector para almacenar mediciones pasadas, use tres colas de prioridad ordenadas (algo así como un montón mínimo / máximo). Estas tres listas dividen la entrada en tres: elementos mayores que la media, elementos menores que la media y elementos iguales a la media.
2. (Casi) cada vez que agrega un elemento, la media cambia, por lo que debemos repartir. Lo crucial es la naturaleza ordenada de las particiones, lo que significa que en lugar de escanear cada elemento de la lista para repartir, solo necesitamos leer los elementos que estamos moviendo. Mientras que en el peor de los casos esto todavía requeriráO(n) operaciones de movimiento, para muchos casos de uso esto no es así.
3. Usando una contabilidad inteligente, podemos asegurarnos de que la desviación se calcule correctamente en todo momento, al repartir y al agregar nuevos elementos.

Algún código de muestra, en python, está debajo. Tenga en cuenta que solo permite agregar elementos a la lista, no eliminarlos. Esto podría agregarse fácilmente, pero en el momento en que escribí esto no lo necesitaba. En lugar de implementar las colas de prioridad yo mismo, he usado la lista ordenada del excelente paquete blist de Daniel Stutzbach , que usa B + Tree s internamente.

Considere este código licenciado bajo la licencia MIT . No se ha optimizado ni pulido significativamente, pero ha funcionado para mí en el pasado. Nuevas versiones estarán disponibles aquí . Avíseme si tiene alguna pregunta o si encuentra algún error.

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 los síntomas persisten, consulta a tu médico.

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

MAD (x) es solo dos cálculos medios concurrentes, cada uno de los cuales se puede hacer en línea a través del binmedian algoritmo .

Puede encontrar el documento asociado, así como el código C y FORTRAN en línea aquí .

(esto es solo el uso de un truco inteligente además del ingenioso truco de Shabbychef, para ahorrar memoria).

Apéndice:

Existe una gran cantidad de métodos antiguos de múltiples pasos para calcular cuantiles. Un enfoque popular es mantener / actualizar un reservorio de observaciones de tamaño determinista seleccionado al azar de la corriente y calcular cuantiles de forma recursiva (ver esta revisión) en este reservorio. Este enfoque (y el relacionado) son reemplazados por el propuesto anteriormente.

Lo siguiente proporciona una aproximación inexacta, aunque la imprecisión dependerá de la distribución de los datos de entrada. Es un algoritmo en línea, pero solo se aproxima a la desviación absoluta. Se basa en un algoritmo bien conocido para calcular la varianza en línea, descrito por Welford en la década de 1960. Su algoritmo, traducido a R, se ve así:

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

Se realiza de manera muy similar a la función de varianza incorporada de 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 

La modificación del algoritmo para calcular la desviación absoluta simplemente implica una sqrtllamada adicional . Sin embargo, sqrtpresenta imprecisiones que se reflejan en el resultado:

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

Los errores, calculados como arriba, son mucho mayores que para el cálculo de la varianza:

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

Sin embargo, dependiendo de su caso de uso, esta magnitud de error puede ser aceptable.