Algorithme de variance roulante

Question

J'essaie de trouver un algorithme efficace et numériquement stable pour calculer une variance glissante (par exemple, une variance sur une fenêtre glissante de 20 périodes). Je connais l'algorithme Algorithme de Welford qui calcule efficacement la variance courante pour un flux de nombres (il ne nécessite qu'un seul passage), mais je ne suis pas sûr que cela puisse être adapté pour une fenêtre glissante. J'aimerais également que la solution évite les problèmes de précision discutés en haut de l'article. cet article par John D. Cook. Une solution dans n'importe quelle langue fait l'affaire.

Answer 1

J'ai également rencontré ce problème. Il existe d'excellents articles sur le calcul de la variance cumulative courante, comme celui de John Cooke. Calculer avec précision la variance courante et le billet de Digital explorations, Code Python pour le calcul des variances, de la covariance et du coefficient de corrélation d'un échantillon et d'une population . Je n'ai pas pu en trouver qui soient adaptés à une fenêtre roulante.

El Écarts types courants Le post de Subluminal Messages a été essentiel pour faire fonctionner la formule de la fenêtre roulante. Jim utilise la puissance de la somme des différences au carré des valeurs, alors que l'approche de Welford consiste à utiliser la somme des différences au carré de la moyenne. La formule est la suivante :

PSA aujourd'hui = PSA(hier) + (((x aujourd'hui * x aujourd'hui) - x hier)) / n

• x = valeur dans votre série temporelle
• n = nombre de valeurs que vous avez analysées jusqu'à présent.

Mais, pour convertir la formule Power Sum Average en une variété fenêtrée, vous devez modifier la formule comme suit :

PSA aujourd'hui = PSA hier + (((x today * x today) - (x yesterday * x Yesterday) / n

• x = valeur dans votre série temporelle
• n = nombre de valeurs que vous avez analysées jusqu'à présent.

Vous aurez également besoin de la formule de la moyenne mobile simple glissante :

SMA aujourd'hui = SMA hier + ((x aujourd'hui - x aujourd'hui - n) / n

• x = valeur dans votre série temporelle
• n = période utilisée pour votre fenêtre de roulement.

À partir de là, vous pouvez calculer la variance de la population mobile :

Population Var aujourd'hui = (PSA aujourd'hui * n - n * SMA aujourd'hui * SMA aujourd'hui) / n

Ou la variance de l'échantillon mobile :

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

J'ai abordé ce sujet avec un exemple de code Python dans un article de blog il y a quelques années, Variation courante .

J'espère que cela vous aidera.

Remarque : j'ai fourni des liens vers tous les articles du blog et les formules mathématiques en Latex (images) pour cette réponse. Mais, en raison de ma faible réputation (< 10) ; je suis limité à seulement 2 hyperliens et absolument aucune image. Désolé pour cela. J'espère que cela n'enlève rien au contenu.

Answer 2

J'ai été confronté au même problème.

La moyenne est simple à calculer de manière itérative, mais vous devez conserver l'historique complet des valeurs dans un tampon circulaire.

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;

J'ai adapté l'algorithme de Welford et il fonctionne pour toutes les valeurs que j'ai testées.

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;

Pour obtenir la variance actuelle, il suffit de diviser varSum par la taille de la fenêtre : variance = varSum / window_size;

Answer 3

Si vous préférez le code aux mots (fortement inspiré du post 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 4

En fait, l'algorithme de Welford peut facilement être adapté pour calculer pondéré Variance. Et en fixant les poids à -1, vous devriez être en mesure d'annuler efficacement les éléments. Je n'ai pas vérifié les mathématiques pour savoir si elles autorisent les poids négatifs, mais à première vue, elles devraient le faire !

J'ai fait une petite expérience en utilisant 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);
}
}

J'obtiens environ ~14 chiffres de précision par rapport à l'algorithme exact à deux passages ; c'est à peu près tout ce que l'on peut attendre des doubles. Notez que Welford fait a un certain coût de calcul en raison des divisions supplémentaires - il prend environ deux fois plus de temps que l'algorithme exact à deux passages. Si la taille de votre fenêtre est petite, il peut être beaucoup plus judicieux de recalculer la moyenne et, dans un second temps, la variance. chaque temps.

J'ai ajouté cette expérience comme test unitaire à ELKI, vous pouvez voir la source complète ici : http://elki.dbs.ifi.lmu.de/browser/elki/trunk/test/de/lmu/ifi/dbs/elki/math/TestSlidingVariance.java il se compare également à la variance exacte à deux passages.

Cependant, sur des ensembles de données asymétriques, le comportement peut être différent. Cet ensemble de données est évidemment uniformément distribué, mais j'ai également essayé un tableau trié et cela a fonctionné.

Mise à jour : nous avons publié un article avec des détails sur différents schémas de pondération pour la (co-)variance :

Schubert, Erich, et Michael Gertz. " Calcul parallèle numériquement stable de la (co-)variance. " Actes de la 30e conférence internationale sur la gestion des bases de données scientifiques et statistiques. ACM, 2018. (A remporté le prix du meilleur article de la SSDBM).

Il est également question de la façon dont la pondération peut être utilisée pour paralléliser le calcul, par exemple avec AVX, les GPU ou sur des grappes.

Answer 5

Voici une approche "diviser pour mieux régner" qui a O(log k) -les mises à jour temporelles, où k est le nombre d'échantillons. Cela devrait être relativement stable pour les mêmes raisons que la sommation par paire et les FFT sont stables, mais c'est un peu compliqué et la constante n'est pas grande.

Supposons que nous ayons une séquence A de longueur m avec une moyenne E(A) et la variance V(A) et une séquence B de longueur n avec une moyenne E(B) et la variance V(B) . Soit C soit la concaténation de A y B . Nous avons

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

Maintenant, placez les éléments dans un arbre rouge-noir, où chaque nœud est décoré avec la moyenne et la variance du sous-arbre enraciné à ce nœud. Insérer à droite, supprimer à gauche. (Puisque nous n'accédons qu'aux extrémités, un arbre splay pourrait être O(1) amorti, mais je suppose que l'amorti est un problème pour votre application). Si k est connu au moment de la compilation, vous pouvez probablement dérouler la boucle interne à la manière de FFTW.