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

InformationsquelleAutor Abiel | 2011-02-28

23

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:
PSA aujourd'hui = PSA(hier) + (((x aujourd'hui * x aujourd'hui) - x hier)) /n
- x = valeur de votre temps à la série
- n = nombre de valeurs que vous avez analysé jusqu'à présent.
Mais, pour convertir la Somme de la Puissance Moyenne de la formule fenêtré variété vous avez besoin d'ajuster la formule suivante:
PSA aujourd'hui = PSA hier + (((x aujourd'hui * x aujourd'hui) - (x-hier * x Hier) /n
- x = valeur de votre temps à la série
- n = nombre de valeurs que vous avez analysé jusqu'à présent.
Vous aurez également besoin de Rouler Moyenne mobile Simple formule:
SMA aujourd'hui = SMA hier + ((x aujourd'hui - x d'aujourd'hui - n) /n
- x = valeur de votre temps à la série
- n = période utilisée pour votre roulant fenêtre.
À partir de là, vous pouvez calculer les Rolling Variance de Population:

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

Ou les Rolling Variance de l'Échantillon:

Échantillon Var aujourdhui = (PSA aujourd'hui * n - n * SMA aujourd'hui * SMA aujourd'hui) /(n - 1)

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.

Veuillez noter: j'ai fourni des liens vers tous les articles de blog et les formules mathématiques
en Latex (images) pour répondre à cette question. Mais, à cause de mon faible réputation (<
10); je suis limité à seulement 2 liens hypertextes et absolument pas d'images. Désolé
à propos de cette. Espérons que ce n'est pas le contenu.
- Dans cette formule: Population Var today = (PSA today * n - n * SMA today * SMA today) / n - pourquoi ne pas supprimer n? Population Var today = (PSA today - SMA today * SMA today).
- En raison de la quadrature des échantillons dans la formule, cet algorithme présente le très numérique inexactitude que l'opération a été d'essayer d'éviter.
- Yup, pas numériquement stable approche. La chose la plus proche à une réponse correcte est par @DanS ci-dessous.
- Merci pour l'explication, voici un C# de mise en œuvre de gist.github.com/mattdot/d459b1cb15480fefd953841a1ac70be8
InformationsquelleAutor Mike Taylor
18

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.
```
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é Welford de l'algorithme, et il fonctionne pour toutes les valeurs que j'ai testé 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 le courant de la variance juste diviser varSum par la taille de la fenêtre: variance = varSum /window_size;
- Il pourrait être un peu plus stable pour varSum += (x_new + x_old - mean - new_mean) * (x_new - x_old), où x_old = xs[next_index], que vous supprimez un large mean * new_mean en additionnant les deux articles que vous soustraire à la mise à jour varSum. 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.
- Afin de clarifier Jaime réponse, il a fait quelques algèbre de prendre DanS la varSum équation et la distribution de la multiplication. certains termes annuler, mais vous aurez également à effectuer l'astuce de mettre en x_new * x_old - x_new * x_old pour arriver à son résultat
- Très tard commentaire: Pourquoi êtes-vous de la plongée par window_size et pas window_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.
InformationsquelleAutor DanS

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

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

InformationsquelleAutor Joachim

5

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 longueur m avec une moyenne de E(A) et de la variance V(A), et une séquence B de longueur n avec une moyenne de E(B) et de la variance V(B). Laissez C être 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, 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.) Si k est connu à la compilation, vous pourriez probablement dérouler la boucle interne FFTW de style.
- (Note: c'est bien de calculer q = 1 - p, sauf si k est stupendously grande.)
- Bon d'accord, c'est essentiellement Chan et coll., s'algorithme parallèle comme décrit sur Wikipédia. C'est ce que je peux avoir pour pas faire défiler vers le bas...
- Pouvez-vous nous expliquer un peu plus en détail comment vous le feriez appliquer cet algorithme à l'écart sur une fenêtre mobile? Je suis un peu familier avec le Chan et coll approche, mais pensé comme un passe-méthode pour le calcul d'une seule variation par rapport à l'ensemble de l'échantillon, avec l'avantage que le problème peut être divisé en parties, qui sont exécutés en parallèle.
- Chan et al a donné un moyen pour calculer les statistiques d'un enchaînement de pièces étant donné les statistiques des parties. Le niveau élevé de l'idée est de maintenir une collection de pièces (en fait juste leurs statistiques) de telle sorte que n'importe quelle fenêtre est la concaténation de O(log k) les pièces. Une façon est avec un arbre binaire équilibré, mais comme Rex points, c'est exagéré, et il nous suffit de maintenir des statistiques pour la aligné pièces dont les tailles sont des puissances de deux (p. ex., [0, 1), [1, 2), [0, 2), [2, 3), [3, 4), [2, 4), [0, 4), etc.)
InformationsquelleAutor userOVER9000
4

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:
```
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);
}
}
```
- 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é.

InformationsquelleAutor Erich Schubert

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.

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

merci pour l'idée de la synthèse en python. Je n'aime pas comment la taille des fenêtres devenir WIN_SIZE - 1dans le cas où le bloc else est entré. Donc, si WIN_SIZE était de 10 lorsque push est appelé et nous ajouter, c'est encore de 10 à cause de la deque constructeur option est utilisée, puis dans le else bloc popleft réduit la taille de plus de 9. Alors peut-être maxlen=WIN_SIZE + 1? Ou ne pas utiliser le maxlen option. Aussi, peut déposer l' n variable et utiliser len(self.windows).
dans le get_var méthode le dénominateur devrait être self.n ou len(self.windows)

InformationsquelleAutor ewerlopes

1

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.

InformationsquelleAutor Andrew White
1

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:
```
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
```
- Je crois que cette solution est sensible aux problèmes de stabilité mentionné dans le lien dans mon premier post (johndcook.com/standard_deviation.html). En particulier, lorsque les valeurs d'entrée et de grandes et leur différence est faible, le résultat pourrait être négatif. Je n'ai aucun contrôle sur l'entrée, donc je préfère éviter cette approche.
- Oh, je vois. Pouvez-vous nous dire à propos de l'entrée? Utilisation prévue? Est-ce un problème que vous pouvez simplement jeter plus de bits à 64 bits float, arbitraire précision arithmétique, etc.)? Les erreurs d'arrondi aller loin si vous trump de la saisie de chiffres significatifs, non?
- d'accord-ce qui a des problèmes de stabilité. Imaginez 1000 échantillons de près de 1 000 000.0, 20 échantillons proche de zéro.
- S: Le roulement de la variance est ce que c'est. Il y a peut être beaucoup de choses dans la transition de 1 million de dollars pour ~zéro, mais c'est la nature de la bête. Que, et la première 980 1000 ~1 million de valeurs de l'image lorsque le changement se produit de toute façon. Mon commentaire a suggéré que, si vous avez eu assez de chiffres significatifs dans vos calculs, rien de tout cela n'aurait de l'importance.
- D'entrée pourrait vraiment être quelque chose. Grandeur de valeur pourrait certainement être dans les milliers de milliards, et alors que les données d'origine aura seulement une précision à quelques décimales, les utilisateurs seront en mesure de transformer leurs données (par exemple en divisant par un scalaire) avant le calcul de la variance.
InformationsquelleAutor John
1

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:
```
|||||||||||||||||||||||||  //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 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.
- "utilisez votre standard E(x^2)-E(x)^2 formule" Non, ne l'a pas; il n'est pas même à distance stable. Adapter l'un des meilleurs algorithmes.
- Pourquoi êtes-vous inquiet au sujet de la stabilité de plus de 20 articles? Cumulatif des erreurs qui s'accumulent au fil des millions d'entrées sont un problème (en particulier lors de la prise de roulement de la fenêtre), mais ce n'est pas la question ici.
- Je suis inquiet parce que c'est un problème. Allez lire l'article de Wikipedia, et si vous n'êtes toujours pas convaincu, essayez le calcul de la variance de 20 iid échantillons de N(1, 1e-10).
- Je n'ai pas vu ce réellement être un problème pour tout réaliste jeu de données avec des unités et de l'origine, mais juste assez, si c'est ce que l'OP veut...
InformationsquelleAutor Rex Kerr
1

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:
- lorsqu'un nouvel élément arrive, supprimez la RunningStat exemple, ajouter l'élément à tous les 20 instances, et incrémenter le "compteur" (modulo 20), qui identifie le nouveau "plein" RunningStat instance
- à tout moment, vous pouvez consulter la "complète" de l'instance pour obtenir votre variante.
Vous 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 dure Mk et Sk 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.

InformationsquelleAutor Matthieu M.
0

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).
```
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 la plus ancienne de l'échantillon dans la fenêtre que vous souhaitez supprimer.

InformationsquelleAutor vibe

Vous devez vous connecter pour publier un commentaire.