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

16

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.

gianluca
fuente

Respuestas:

6

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 (MM=1000Bi,ji=1Mj=1NNkBi,jj). 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 ](i1)/Mi/MN[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 = 1Kjy elelementojde la desviación absoluta media se estima en1

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

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.

shabbychef
fuente
Wow, enfoque realmente interesante. No sabía sobre eso, y lo tendré en cuenta. Desafortunadamente, en este caso no funcionará, ya que tengo requisitos realmente restrictivos desde el punto de vista del uso de memoria, por lo que M debería ser realmente pequeño y supongo que habría demasiada pérdida de precisión.
gianluca
@gianluca: parece que tiene 1. muchos datos, 2. recursos de memoria limitados, 3. requisitos de alta precisión. ¡Puedo ver por qué este problema te está volviendo loco! Tal vez, como lo menciona @kwak, puede calcular alguna otra medida de propagación: MAD, IQR, desviación estándar. Todos ellos tienen enfoques que podrían funcionar para su problema.
shabbychef
gianluca:> Danos una idea más cuantitativa sobre el tamaño de la memoria, las matrices y la precisión que deseas. Sin embargo, es muy posible que su pregunta se responda mejor en stackoverflow.
usuario603
4

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.

fmark
fuente
2
Me falta algo: si tiene que "mantener todo el vector en la memoria", ¿cómo califica esto como un algoritmo "en línea"?
whuber
@whuber No, no falta algo, supongo que no es un algoritmo en línea. Requiere O(n)memoria y, en el peor de los casos, toma O (n) tiempo para cada elemento agregado. Sin embargo, en los datos distribuidos normalmente (y probablemente en otras distribuciones) funciona de manera bastante eficiente.
fmark el
3

XXXss2/π

shabbychef
fuente
Esa es una idea interesante. Tal vez podría complementarlo con la detección en línea de valores atípicos y utilizarlos para modificar la estimación a medida que avanza.
whuber
Probablemente podría usar el método de Welford para calcular la desviación estándar en línea que documenté en mi segunda respuesta.
fmark el
1
Sin embargo, se debe tener en cuenta que de esta manera se podría perder la robustez de los estimadores, como el MAD explícito, que a veces conducen a su elección frente a alternativas más simples.
Quartz
2

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.

usuario603
fuente
¿Podría detallar o hacer referencia a la relación entre MAD y las dos medianas?
Quartz
medi=1n|ximedi=1n|
Hehm, en realidad quise decir si puedes explicar cómo esta relación permite que las dos medianas sean concurrentes; me parecen dependientes, ya que las entradas a la mediana externa pueden cambiar en cada muestra agregada al cálculo interno. ¿Cómo los realizarías en paralelo?
Quartz
Tengo que volver al papel binmedian para obtener más detalles ... pero dado un valor calculado de la mediana (metromireyo=1norteXyo) y un nuevo valor de Xnorte+1 el algoritmo podría calcular metromireyo=1norte+1Xyo mucho más rápido que O(norte) identificando el contenedor al que Xnorte+1pertenece No veo cómo esta idea no podría generalizarse a la mediana externa en el cómputo loco.
user603
1

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.

historgram of differences

fmark
fuente
Esto no da la respuesta exacta, por la siguiente razón: yoXyoyoXyo. Estás calculando lo primero, mientras que el OP quiere lo segundo.
shabbychef
Estoy de acuerdo en que el método es inexacto. Sin embargo, no estoy de acuerdo con su diagnóstico de la inexactitud. El método de Welford para calcular la varianza, que ni siquiera contiene un sqrt, tiene un error similar. Sin embargo, a medida que se nhace más grande, error/nse desvanece muy poco, sorprendentemente rápido.
fmark
Welford's method has no sqrt because it is computing the variance, not the standard deviation. By taking the sqrt, it seems like you are estimating the standard deviation, not the mean absolute deviation. am I missing something?
shabbychef
@shabbychef Each iteration of Welfords is calculating the contribution of the new datapoint to the absolute deviation, squared. So I take the square root of each contribution squared to get back to the absolute deviance. You might note, for example, that I take the square root of the delta before I add it to the deviance sum, rather than afterward as in the case of the standard deviation.
fmark
3
I see the problem; Welfords obscures the problem with this method: the online estimate of the mean is being used instead of the final estimate of the mean. While Welford's method is exact (up to roundoff) for variance, this method is not. The problem is not due to the sqrt imprecision. It is because it uses the running mean estimate. To see when this will break, try xs <- sort(rnorm(n.testitems)) When I try this with your code (after fixing it to return a.dev / n), I get relative errors on the order of 9%-16%. So this method is not permutation invariant, which could cause havoc...
shabbychef