algorithm - keywords - meta tags seo 2018

Algoritmo de varianza rodante

Question

Estoy tratando de encontrar un algoritmo eficiente, numéricamente estable para calcular una varianza continua (por ejemplo, una varianza en una ventana rotativa de 20 periodos). Conozco el algoritmo de Welford que calcula eficientemente la varianza de ejecución para una secuencia de números (solo requiere una pasada), pero no estoy seguro si esto se puede adaptar para una ventana móvil. También me gustaría la solución para evitar los problemas de precisión discutidos en la parte superior de este artículo por John D. Cook. Una solución en cualquier idioma está bien.

Answer 1

Aquí hay otra solución O(log k) : encuentra cuadrados la secuencia original, luego suma pares, luego cuadruplica, etc. (Necesitarás un poco de amortiguación para poder encontrar todos estos eficientemente). Luego sume esos valores que necesita para obtener su respuesta. Por ejemplo:

|||||||||||||||||||||||||  // Squares
| | | | | | | | | | | | |  // Sum of squares for pairs
|   |   |   |   |   |   |  // Pairs of pairs
|       |       |       |  // (etc.)
|               |
   ^------------------^    // Want these 20, which you can get with
        |       |          // one...
    |   |       |   |      // two, three...
                    | |    // four...
   ||                      // five stored values.

Ahora usa su fórmula estándar E (x ^ 2) -E (x) ^ 2 y listo. (No, si necesita una buena estabilidad para conjuntos pequeños de números; esto suponía que era solo la acumulación de error de desplazamiento lo que causaba problemas).

Dicho eso, sumar 20 números al cuadrado es muy rápido actualmente en la mayoría de las arquitecturas. Si estuvieras haciendo más, digamos, un par de cientos, un método más eficiente sería claramente mejor. Pero no estoy seguro de que la fuerza bruta no sea el camino a seguir aquí.

Answer 2

Aquí hay un enfoque de dividir y conquistar que tiene actualizaciones O(log k) -time, donde k es el número de muestras. Debería ser relativamente estable por las mismas razones que la suma por pares y las FFT son estables, pero es un poco complicado y la constante no es buena.

Supongamos que tenemos una secuencia A de longitud m con media E(A) y varianza V(A) , y una secuencia B de longitud n con media E(B) y varianza V(B) . Deje C ser la concatenación de A y B Tenemos

p = m / (m + n)
q = n / (m + n)
E(C) = p * E(A) + q * E(B)
V(C) = p * (V(A) + (E(A) + E(C)) * (E(A) - E(C))) + q * (V(B) + (E(B) + E(C)) * (E(B) - E(C)))

Ahora, rellena los elementos en un árbol rojo-negro, donde cada nodo está decorado con la media y la varianza del subárbol enraizado en ese nodo. Insertar a la derecha; eliminar a la izquierda. (Dado que solo estamos accediendo a los extremos, un árbol de distribución podría ser O(1) amortizado, pero supongo que amortizado es un problema para su aplicación). Si se conoce k en tiempo de compilación, probablemente podría desenrollar el interior loop estilo FFTW.

Answer 3

En realidad, el algoritmo Welfords puede adaptarse fácilmente para calcular la varianza ponderada . Y al establecer los pesos en -1, debería poder cancelar los elementos de manera efectiva. No he comprobado las matemáticas si permite pesos negativos, ¡pero a primera vista debería!

Realicé un pequeño experimento usando ELKI :

void testSlidingWindowVariance() {
MeanVariance mv = new MeanVariance(); // ELKI implementation of weighted Welford!
MeanVariance mc = new MeanVariance(); // Control.

Random r = new Random();
double[] data = new double[1000];
for (int i = 0; i < data.length; i++) {
  data[i] = r.nextDouble();
}

// Pre-roll:
for (int i = 0; i < 10; i++) {
  mv.put(data[i]);
}
// Compare to window approach
for (int i = 10; i < data.length; i++) {
  mv.put(data[i-10], -1.); // Remove
  mv.put(data[i]);
  mc.reset(); // Reset statistics
  for (int j = i - 9; j <= i; j++) {
    mc.put(data[j]);
  }
  assertEquals("Variance does not agree.", mv.getSampleVariance(),
    mc.getSampleVariance(), 1e-14);
}
}

Obtengo ~ 14 dígitos de precisión en comparación con el algoritmo exacto de dos pasadas; esto es casi todo lo que se puede esperar de los dobles. Tenga en cuenta que Welford sí tiene algún costo computacional debido a las divisiones adicionales; tarda aproximadamente el doble que el algoritmo exacto de dos pasadas. Si el tamaño de su ventana es pequeño, puede ser mucho más sensato recalcular la media y luego, en un segundo pase, la varianza cada vez.

He añadido este experimento como prueba unitaria a ELKI, puedes ver la fuente completa aquí: http://elki.dbs.ifi.lmu.de/browser/elki/trunk/test/de/lmu/ifi/dbs/elki/math/TestSlidingVariance.java también se compara con la varianza de dos pasadas exacta.

Sin embargo, en conjuntos de datos asimétricos, el comportamiento puede ser diferente. Este conjunto de datos obviamente está distribuido uniformemente; pero también probé una matriz ordenada y funcionó.

Answer 4

Espero que se demuestre lo contrario, pero no creo que esto pueda hacerse "rápidamente". Dicho eso, una gran parte del cálculo es mantener un registro del EV sobre la ventana que se puede hacer fácilmente.

Me iré con la pregunta: ¿estás seguro de que necesitas una función de ventana? A menos que trabaje con ventanas muy grandes, probablemente sea mejor usar un algoritmo predefinido bien conocido.

Answer 5

Esto es solo una pequeña adición a la excelente respuesta proporcionada por DanS. Las siguientes ecuaciones son para eliminar la muestra más antigua de la ventana y actualizar la media y la varianza. Esto es útil, por ejemplo, si desea tomar ventanas más pequeñas cerca del borde derecho de su flujo de datos de entrada (es decir, simplemente elimine la muestra de la ventana más antigua sin agregar una nueva muestra).

window_size -= 1; % decrease window size by 1 sample
new_mean = prev_mean + (prev_mean - x_old) / window_size
varSum = varSum - (prev_mean - x_old) * (new_mean - x_old)

Aquí, x_old es la muestra más antigua en la ventana que desea eliminar.

Answer 6

He estado lidiando con el mismo problema.

La media es simple de calcular iterativamente, pero necesita mantener el historial completo de valores en un buffer circular.

next_index = (index + 1) % window_size;    // oldest x value is at next_index, wrapping if necessary.

new_mean = mean + (x_new - xs[next_index])/window_size;

He adaptado el algoritmo de Welford y funciona para todos los valores con los que he probado.

varSum = var_sum + (x_new - mean) * (x_new - new_mean) - (xs[next_index] - mean) * (xs[next_index] - new_mean);

xs[next_index] = x_new;
index = next_index;

Para obtener la varianza actual simplemente divida varSum por el tamaño de la ventana: variance = varSum / window_size;

Answer 7

Me he topado con este problema también. Hay algunas publicaciones excelentes en el cálculo de la varianza acumulada en ejecución, como la publicación de la varianza en ejecución precisa de John Cooke y la publicación de exploraciones digitales, el código de Python para calcular varianzas de muestra y población, covarianza y coeficiente de correlación . Simplemente no pude encontrar ninguno que se haya adaptado a una ventana rodante.

La publicación Desviaciones estándar en curso por mensajes subluminales fue fundamental para que la fórmula de la ventana móvil funcionara. Jim toma la suma de poder de las diferencias al cuadrado de los valores versus el enfoque de Welford de usar la suma de las diferencias al cuadrado de la media. Fórmula de la siguiente manera:

PSA hoy = PSA (ayer) + (((x hoy * x hoy) - x ayer)) / n

• x = valor en tu serie temporal
• n = número de valores que ha analizado hasta ahora.

Pero, para convertir la fórmula del promedio de la suma de energía a una variedad con ventana, debe ajustar la fórmula a lo siguiente:

PSA hoy = PSA ayer + (((x hoy * x hoy) - (x ayer * x ayer) / n

• x = valor en tu serie temporal
• n = número de valores que ha analizado hasta ahora.

También necesitará la fórmula Rolling Simple Moving Average:

SMA hoy = SMA ayer + ((x hoy - x hoy - n) / n

• x = valor en tu serie temporal
• n = período utilizado para su ventana móvil.

A partir de ahí, puede calcular la variación de la población móvil:

Población Var hoy = (PSA hoy * n - n * SMA hoy * SMA hoy) / n

O la varianza de la muestra rodante:

Sample Var today = (PSA today * n - n * SMA today * SMA today) / (n - 1)

He cubierto este tema junto con el código de Python de ejemplo en una publicación de blog hace unos años, Running Variance .

Espero que esto ayude.

Tenga en cuenta: proporcioné enlaces a todas las publicaciones de blog y fórmulas matemáticas en Latex (imágenes) para esta respuesta. Pero, debido a mi baja reputación (<10); Estoy limitado a solo 2 hipervínculos y absolutamente ninguna imagen. Perdón por esto. Espero que esto no quite el contenido.

Answer 8

Por solo 20 valores, es trivial adaptar el método expuesto aquí (aunque no dije rápido).

Simplemente puede elegir una matriz de 20 de estas clases RunningStat .

Los primeros 20 elementos de la transmisión son algo especiales, sin embargo, una vez hecho esto, es mucho más simple:

• cuando llega un nuevo elemento, borre la instancia actual de RunningStat , agregue el elemento a las 20 instancias e incremente el "contador" (módulo 20) que identifica la nueva instancia "completa" de RunningStat
• en cualquier momento, puede consultar la instancia actual "completa" para obtener su variante de ejecución.

Obviamente notará que este enfoque no es realmente escalable ...

También puede observar que hay cierta redudancia en los números que mantenemos (si va con la clase completa de RunningStat ). Una mejora obvia sería mantener las 20 duraciones Mk y Sk directamente.

No puedo pensar en una fórmula mejor usando este algoritmo particular, me temo que su formulación recursiva de alguna manera ata nuestras manos.

Answer 9

Sé que esta pregunta es antigua, pero en caso de que alguien más esté interesado aquí sigue el código python. Está inspirado en johndcook blog post, @ Joachim''s, @ DanS''s code y @Jaime comments. El siguiente código todavía da pequeñas imprecisiones para tamaños de ventanas de datos pequeños. Disfrutar.

from __future__ import division
import collections
import math


class RunningStats:
    def __init__(self, WIN_SIZE=20):
        self.n = 0
        self.mean = 0
        self.run_var = 0
        self.WIN_SIZE = WIN_SIZE

        self.windows = collections.deque(maxlen=WIN_SIZE)

    def clear(self):
        self.n = 0
        self.windows.clear()

    def push(self, x):

        self.windows.append(x)

        if self.n <= self.WIN_SIZE:
            # Calculating first variance
            self.n += 1
            delta = x - self.mean
            self.mean += delta / self.n
            self.run_var += delta * (x - self.mean)
        else:
            # Adjusting variance
            x_removed = self.windows.popleft()
            old_m = self.mean
            self.mean += (x - x_removed) / self.WIN_SIZE
            self.run_var += (x + x_removed - old_m - self.mean) * (x - x_removed)

    def get_mean(self):
        return self.mean if self.n else 0.0

    def get_var(self):
        return self.run_var / (self.WIN_SIZE - 1) if self.n > 1 else 0.0

    def get_std(self):
        return math.sqrt(self.get_var())

    def get_all(self):
        return list(self.windows)

    def __str__(self):
        return "Current window values: {}".format(list(self.windows))

Answer 10

Si prefiere el código sobre las palabras (basado en gran medida en la publicación de DanS): http://calcandstuff.blogspot.se/2014/02/rolling-variance-calculation.html

public IEnumerable RollingSampleVariance(IEnumerable data, int sampleSize)
{
    double mean = 0;
    double accVar = 0;

    int n = 0;
    var queue = new Queue(sampleSize);

    foreach(var observation in data)
    {
        queue.Enqueue(observation);
        if (n < sampleSize)
        {
            // Calculating first variance
            n++;
            double delta = observation - mean;
            mean += delta / n;
            accVar += delta * (observation - mean);
        }
        else
        {
            // Adjusting variance
            double then = queue.Dequeue();
            double prevMean = mean;
            mean += (observation - then) / sampleSize;
            accVar += (observation - prevMean) * (observation - mean) - (then - prevMean) * (then - mean);
        }

        if (n == sampleSize)
            yield return accVar / (sampleSize - 1);
    }
}

Answer 11

Supongo que hacer un seguimiento de tus 20 muestras, Suma (X ^ 2 de 1..20), y Suma (X de 1..20) y luego volver a calcular sucesivamente las dos sumas en cada iteración no es lo suficientemente eficiente. Es posible recalcular la nueva varianza sin agregar, cuadrar, etc., todas las muestras cada vez.

Como en:

Sum(X^2 from 2..21) = Sum(X^2 from 1..20) - X_1^2 + X_21^2
Sum(X from 2..21) = Sum(X from 1..20) - X_1 + X_21