Algorithme de variance de roulement

J’essaie de trouver un algorithme efficace et numériquement stable pour calculer une variance de roulement (par exemple, une variance sur une fenêtre glissante de 20 périodes). Je connais l’ algorithme de Welford qui calcule efficacement la variance en cours d’exécution pour un stream de nombres (il ne nécessite qu’un seul passage), mais je ne suis pas sûr que cela puisse être adapté à une fenêtre en cours. J’aimerais aussi que la solution évite les problèmes d’exactitude discutés en haut de cet article de John D. Cook. Une solution dans n’importe quelle langue convient.

J’ai également rencontré ce problème. Il existe d’excellents articles sur le calcul de la variance cumulative courante tels que le calcul de la variance en cours avec calcul précis de John Cooke et le post des explorations numériques, le code Python pour calculer les variances de l’échantillon et de la population, la covariance et le coefficient de corrélation . Juste ne pouvait pas en trouver qui ont été adaptés à une fenêtre roulante.

Le post sur les écarts types en cours d’exécution par messages sublétaux était essentiel pour que la formule de la fenêtre mobile fonctionne. Jim prend la sum des différences au carré des valeurs par rapport à la méthode de Welford consistant à utiliser la sum des différences au carré de la moyenne. Formule comme suit:

PSA aujourd’hui = PSA (hier) + (((x aujourd’hui * x aujourd’hui) – x hier)) / n

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

Cependant, 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 aujourd’hui * x aujourd’hui) – (x hier * x hier) / n

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

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

SMA aujourd’hui = SMA hier + ((x aujourd’hui – x aujourd’hui – n) / n

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

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 roulant:

Sample Var today = (PSA aujourd’hui * n – n * SMA aujourd’hui * SMA aujourd’hui) / (n – 1)

J’ai abordé ce sujet avec un exemple de code Python dans un article de blog il y a quelques années, Running Variance .

J’espère que cela t’aides.

S’il vous plaît noter: J’ai fourni des liens vers tous les billets de 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 ça. J'espère que cela n'enlève rien au contenu.

Je me suis occupé du même problème.

Mean est simple à calculer itérativement, 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 cela fonctionne pour toutes les valeurs que j’ai testées avec.

 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, divisez simplement varSum par la taille de la fenêtre: variance = varSum / window_size;

Si vous préférez le code sur les mots (fortement basé sur le 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); } } 

Voici une approche de diviser et conquérir qui a des mises à jour O(log k) -time, où k est le nombre d’échantillons. Il devrait être relativement stable pour les mêmes raisons que la sum des paires et les FFT sont stables, mais c’est un peu compliqué et la constante n’est pas géniale.

Supposons que nous ayons une séquence A de longueur m avec la moyenne E(A) et la variance V(A) et une séquence B de longueur n avec la moyenne E(B) et la variance V(B) . Soit C la concaténation de A et 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, remplissez 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. (Comme nous accédons uniquement aux extrémités, un arbre splay peut être O(1) amorti, mais je suppose que l’amortissement est un problème pour votre application.) Si k est connu à la compilation, vous pourriez probablement dérouler le fichier interne boucle de style FFTW.

En fait, l’algorithme de Welfords peut facilement être adapté à l’AFAICT pour calculer la variance pondérée . Et en définissant des poids à -1, vous devriez pouvoir annuler efficacement des éléments. Je n’ai pas vérifié si cela permettait des poids négatifs, mais au premier abord, ça devrait l’être!

J’ai effectué une petite expérience avec 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 passes; c'est à peu près autant qu'on peut s'y attendre des doubles. Notez que Welford a un coût de calcul en raison des divisions supplémentaires - cela prend environ deux fois plus de temps que l'algorithme à deux passes exact. Si la taille de votre fenêtre est petite, il peut être beaucoup plus judicieux de recalculer la moyenne et ensuite, dans un second temps, d'écarter la variance.

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 compare également à la variance exacte en deux passes.

Cependant, sur des ensembles de données asymésortingques, le comportement peut être différent. Cet dataset est évidemment dissortingbué uniformément; mais j'ai aussi essayé un tableau sortingé et ça a fonctionné.

Je sais que cette question est ancienne, mais si quelqu’un d’autre est intéressé, suivez le code python. Il s’inspire des articles de blog de johndcook , de @ Joachim, du code de @ DanS et des commentaires de @Jaime. Le code ci-dessous donne toujours de petites imprécisions pour les petites tailles de fenêtres de données. Prendre plaisir.

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

Je suis impatient de prouver que c’est faux, mais je ne pense pas que cela puisse être fait “rapidement”. Cela dit, une grande partie du calcul consiste à garder une trace du véhicule élecsortingque au-dessus de la fenêtre, ce qui peut être fait facilement.

Je vais partir avec la question: êtes-vous sûr d’ avoir besoin d’ une fonction fenêtrée? Sauf si vous travaillez avec de très grandes fenêtres, il est probablement préférable d’utiliser un algorithme prédéfini bien connu.

Je suppose que garder la trace de vos 20 échantillons, Sum (X ^ 2 de 1..20), et Sum (X from 1..20) et ensuite recalculer successivement les deux sums à chaque itération n’est pas assez efficace? Il est possible de recalculer la nouvelle variance sans additionner, mettre au carré, etc., tous les échantillons à chaque fois.

Un péché:

 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 

Voici une autre solution O(log k) : trouvez les carrés de la séquence d’origine, puis additionnez les paires, puis les quadruples, etc. (vous aurez besoin d’un peu de tampon pour pouvoir les retrouver tous efficacement). les valeurs dont vous avez besoin pour obtenir votre réponse. Par exemple:

 

| // 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.

Maintenant, vous utilisez votre formule standard E (x ^ 2) -E (x) ^ 2 et vous avez terminé. (Pas si vous avez besoin d’une bonne stabilité pour de petits nombres de nombres, cela supposait que ce n’était que l’accumulation des erreurs tournantes qui causait des problèmes.)

Cela dit, la sum de 20 nombres au carré est très rapide ces jours-ci sur la plupart des architectures. Si vous faisiez plus, disons quelques centaines, une méthode plus efficace serait clairement meilleure. Mais je ne suis pas sûr que la force brutale ne soit pas la voie à suivre.

Pour seulement 20 valeurs, il est sortingvial d’adapter la méthode exposée ici (je n’ai cependant pas dit vite).

Vous pouvez simplement choisir un tableau de 20 de ces classes RunningStat .

Les 20 premiers éléments du stream sont quelque peu spéciaux, mais une fois cela fait, c’est beaucoup plus simple:

  • lorsqu’un nouvel élément arrive, effacez l’instance RunningStat actuelle, ajoutez l’élément à toutes les 20 instances et incrémentez le “compteur” (modulo 20) qui identifie la nouvelle instance RunningStat “complète”
  • à tout moment, vous pouvez consulter l’instance actuelle “complète” pour obtenir votre variante en cours d’exécution.

Vous remarquerez évidemment que cette approche n’est pas vraiment évolutive …

Vous pouvez également noter qu’il y a une certaine redondance dans les chiffres que nous conservons (si vous RunningStat la classe complète RunningStat ). Une amélioration évidente serait de garder les 20 derniers Mk et Sk directement.

Je ne peux pas penser à une meilleure formule utilisant cet algorithme particulier, je crains que sa formulation récursive ne lie nos mains.

Ceci est juste un ajout mineur à l’excellente réponse fournie par DanS. Les équations suivantes permettent de retirer l’échantillon le plus ancien de la fenêtre et de mettre à jour la moyenne et la variance. Ceci est utile, par exemple, si vous voulez prendre des fenêtres plus petites près du bord droit de votre stream de données en entrée (c.-à-d. Simplement supprimer le plus ancien exemple de fenêtre sans append un nouvel échantillon).

 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) 

Ici, x_old est l’échantillon le plus ancien de la fenêtre que vous souhaitez supprimer.