Rouler la variance de l'algorithme
J'essaie de trouver une façon efficace, numériquement stable algorithme pour calculer un roulement de variance (par exemple, un écart de plus de 20 période roulant fenêtre). Je suis conscient de la Welford algorithme efficacement calcule l'exécution de la variance pour un flux de nombres (elle ne nécessite qu'un seul passage), mais je ne suis pas sûr si cela peut être adapté pour un mobile de la fenêtre. Je voudrais également que la solution pour éviter les problèmes de précision discuté au sommet de cet article par John D. Cook. Une solution dans n'importe quelle langue est fine.
- +1 pour la mention de Welford algorithme; je savais que c'était dans Knuth, mais n'a jamais su la source d'origine
- Bonjour, qu'avez-vous de faire? Avez-vous adapter Chan de l'algorithme? Btw, ne devrait pas kahan somme être en mesure de surmonter les instabilités numériques lors de l'utilisation de la "naïve" de l'approche (suivi la somme des valeurs, et leurs carrés)?
Vous devez vous connecter pour publier un commentaire.
J'ai couru à travers ce type de problème. Il ya quelques grands postes dans le calcul de l'exécution cumulé de la variance comme John Cooke avec Précision le calcul de l'exécution de la variance post et le post de la part de Digital explorations, code Python pour le calcul de l'échantillon et de la population, les écarts, la covariance et le coefficient de corrélation. Seulement ne pouvait pas trouver un qui ont été adaptées à un mobile de la fenêtre.
La L'Exécution Des Écarts-Types post par Subluminal Messages a été critique dans l'obtention de l'roulant fenêtre de la formule de travail. Jim prend le pouvoir somme des carrés des différences des valeurs de rapport de Welford l'approche de l'aide de la somme des différences au carré de la moyenne. Formule comme suit:
Mais, pour convertir la Somme de la Puissance Moyenne de la formule fenêtré variété vous avez besoin d'ajuster la formule suivante:
Vous aurez également besoin de Rouler Moyenne mobile Simple formule:
À partir de là, vous pouvez calculer les Rolling Variance de Population:
Ou les Rolling Variance de l'Échantillon:
J'ai abordé ce sujet ainsi que des exemples de code Python dans un billet de blog quelques années en arrière, L'Exécution De La Variance.
Espère que cette aide.
Population Var today = (PSA today * n - n * SMA today * SMA today) / n
- pourquoi ne pas supprimern
?Population Var today = (PSA today - SMA today * SMA today)
.J'ai eu affaire avec le même problème.
Veux dire, c'est simple à calculer de manière itérative, mais vous devez garder à l'historique complet de valeurs dans une mémoire tampon circulaire.
J'ai adapté Welford de l'algorithme, et il fonctionne pour toutes les valeurs que j'ai testé avec.
Pour obtenir le courant de la variance juste diviser varSum par la taille de la fenêtre:
variance = varSum /window_size;
varSum += (x_new + x_old - mean - new_mean) * (x_new - x_old)
, oùx_old = xs[next_index]
, que vous supprimez un largemean * new_mean
en additionnant les deux articles que vous soustraire à la mise à jourvarSum
. Autre que cela, c'est le plus correct des réponses ici, et c'est dommage que ça n'a pas reçu plus d'amour.varSum
équation et la distribution de la multiplication. certains termes annuler, mais vous aurez également à effectuer l'astuce de mettre enx_new * x_old - x_new * x_old
pour arriver à son résultatwindow_size
et paswindow_size-1
. En d'autres termes: Pourquoi n'êtes-vous pas à l'aide de Bessel de correction. Je remarque que John D. Cook ne comprennent de Bessel de correction dans sa course variance code.Si vous préférez code sur les mots (fortement basé sur les DanS la " post):
http://calcandstuff.blogspot.se/2014/02/rolling-variance-calculation.html
Voici un diviser et conquérir approche qui a
O(log k)
mises à jour en temps, oùk
est le nombre d'échantillons. Il devrait être relativement stable pour les mêmes raisons que les paires de sommation et de la Fft sont stables, mais c'est un peu compliqué et la constante n'est pas très grande.Supposons que nous avons une séquence
A
de longueurm
avec une moyenne deE(A)
et de la varianceV(A)
, et une séquenceB
de longueurn
avec une moyenne deE(B)
et de la varianceV(B)
. LaissezC
être la concaténation deA
etB
. Nous avonsMaintenant, des choses que les éléments en rouge-noir arbre où chaque nœud est décoré avec de la moyenne et de la variance de la sous-arbre enraciné au niveau de ce nœud. Insérer sur le droit; supprimer sur la gauche. (Puisque nous sommes seuls accéder aux extrémités, un splay tree pourrait être
O(1)
amorti, mais je devine amorti est un problème pour votre application.) Sik
est connu à la compilation, vous pourriez probablement dérouler la boucle interne FFTW de style.Fait Welfords algorithme peut AFAICT facilement être adapté pour calculer pondérée de la Variance.
Et par réglage de poids à -1, vous devriez être capable d'annuler les éléments. Je havn'pas vérifié le calcul si le poids est négatif, mais au premier regard il faut!
Je n'ai effectuer une petite expérience à l'aide de ELKI:
- Je obtenir autour de ~14 chiffres de précision par rapport à l'exacte algorithme de deux passes; il s'agit autant que possible de doubles. Notez que Welford ne venir à un certain coût de calcul en raison de la plus-divisions - il faut environ deux fois plus long que l'exacte algorithme de deux passes. Si votre taille de la fenêtre est petite, il peut être beaucoup plus sensible à fait recalculer la moyenne et la puis dans une deuxième passe de la variance chaque temps.
J'ai ajouté cette expérience comme test de l'unité de ELKI, vous pouvez voir l'intégralité du code source ici: http://elki.dbs.ifi.lmu.de/browser/elki/trunk/test/de/lmu/ifi/dbs/elki/math/TestSlidingVariance.java
il compare également à l'exacte deux passes de la variance.
Cependant, inclinées ensembles de données, le comportement peut être différent. Cet ensemble de données est évidemment uniforme distribués; mais j'ai aussi essayé un tableau trié et cela a fonctionné.
Je sais que cette question est ancienne, mais au cas où quelqu'un d'autre est intéressé voici le code python. Il est inspiré par johndcook blog, @Joachim, @DanS le code et @Jaime commentaires. Le code ci-dessous donne toujours des petites imprécisions de données de petite taille, les tailles de fenêtres. Profitez de.
WIN_SIZE - 1
dans le cas où le bloc else est entré. Donc, siWIN_SIZE
était de 10 lorsquepush
est appelé et nous ajouter, c'est encore de 10 à cause de la deque constructeur option est utilisée, puis dans leelse
blocpopleft
réduit la taille de plus de 9. Alors peut-êtremaxlen=WIN_SIZE + 1
? Ou ne pas utiliser lemaxlen
option. Aussi, peut déposer l'n
variable et utiliserlen(self.windows)
.get_var
méthode le dénominateur devrait êtreself.n
oulen(self.windows)
J'ai hâte de le tromper mais je ne pense pas que cela peut être fait "rapidement". Cela dit, une grande partie du calcul est de garder la trace de l'EV-dessus de la fenêtre qui peut être fait facilement.
Je vous laisse avec la question: êtes-vous sûr de besoin fenêtré fonction? Sauf si vous travaillez avec de très grandes fenêtres, il est probablement préférable d'utiliser un bien connu algorithme prédéfini.
Je pense garder une trace de vos 20 échantillons, Sum(X^2 à partir de 1..20) et Sum(X à partir de 1..20) et puis, successivement, à recalculer les deux sommes à chaque itération n'est pas assez efficace? Il est possible de recalculer la nouvelle variance sans l'addition, la quadrature, etc., tous les échantillons à chaque fois.
Comme dans:
Voici un autre
O(log k)
solution: trouver des places, la séquence d'origine, alors la somme des paires, puis quadruple, etc.. (Vous aurez besoin d'un peu d'une mémoire tampon pour être en mesure de trouver toutes ces efficacement.) Puis ajouter ces valeurs que vous avez besoin pour obtenir votre réponse. Par exemple:Maintenant, vous utilisez votre standard E(x^2)-E(x)^2 formule et vous avez terminé.(Pas si vous avez besoin d'une bonne stabilité pour les petits ensembles de nombres; c'est en supposant que c'était seulement de l'accumulation de roulement erreur qui a été à l'origine de problèmes.)Cela dit, résumant 20 carrés des nombres est très rapide ces jours sur la plupart des architectures. Si vous faisiez plus de, disons, quelques centaines de--la méthode la plus efficace serait nettement mieux. Mais je ne suis pas sûr que la force brute n'est pas la voie à suivre ici.
Pour seulement 20 valeurs, il est trivial d'adapter la méthode exposée ici (je n'ai pas dit rapide, cependant).
Vous pouvez tout simplement ramasser un tableau de 20 de ces
RunningStat
classes.Les 20 premiers éléments de la rivière sont un peu spéciale, cependant une fois que c'est fait, c'est beaucoup plus simple:
RunningStat
exemple, ajouter l'élément à tous les 20 instances, et incrémenter le "compteur" (modulo 20), qui identifie le nouveau "plein"RunningStat
instanceVous permettra évidemment de noter que cette approche n'est pas vraiment évolutive...
Vous pouvez également noter qu'il y a quelques redudancy dans les chiffres que nous garder (si vous y allez avec le
RunningStat
complet de la classe). Une amélioration évidente serait de garder le 20 dureMk
etSk
directement.Je ne peux pas penser à une meilleure formule à l'aide de cet algorithme, j'ai peur que sa formulation récursive quelque peu les liens entre nos mains.
C'est juste un mineur de plus de l'excellente réponse fournie par la DanS. Les équations suivantes sont pour la suppression de la plus ancienne de l'échantillon à partir de la fenêtre et la mise à jour de la moyenne et de la variance. Ceci est utile, par exemple, si vous voulez prendre de plus petites fenêtres près du bord droit de votre flux de données d'entrée (c'est à dire il suffit de retirer le plus ancien de la fenêtre de l'échantillon, sans l'ajout d'un nouvel échantillon).
Ici, x_old est la plus ancienne de l'échantillon dans la fenêtre que vous souhaitez supprimer.