Algorithme de calcul de la variance
Les algorithmes de calcul de la variance jouent un rôle majeur dans les statistiques numériques. Une difficulté clé dans la conception de bons algorithmes pour ce problème est que les formules de calcul de la variance peuvent impliquer le calcul de sommes de carrés, qui peuvent provoquer des instabilités numériques de même que des dépassements arithmétiques quand de grandes valeurs apparaissent..
Algorithme naïf
[modifier | modifier le code]Une formule usuelle pour le calcul de la variance d'une population de taille N est :
Avec la correction de Bessel qui permet d'avoir un estimateur sans biais de la variance pour un échantillon fini de n observations, the formula is:
Ainsi, un algorithme naïve algorithm to calculate the estimated variance is given by the following:
- Soient n ← 0, Sum ← 0, SumSq ← 0
- Pour chaque donné x :
- n ← n + 1
- Sum ← Sum + x
- SumSq ← SumSq + x × x
- Var = (SumSq − (Sum × Sum) / n) / (n − 1)
Cet algorithme peut facilement être adapté au calcul de la variance d'une population finie, en divisant par n au lieu de n − 1 à la dernière étape.
Puisque SumSq et (Sum × Sum)/n peuvent être des nombres très proches, une annulation catastrophique peut entrainer une précision du résultat inférieure à la précision machine utilisée pour le calcul. Ainsi, cet algorithme ne doit pas être utilisé en pratique[1],[2], aussi plusieurs algorithmes alternatifs, numériquement stables, ont été proposés[3]. C'est notablement problématique si l'écart type est petit par rapport à la moyenne.
Calcul avec des données décalées
[modifier | modifier le code]La variance est invariante par translation, donc insensible à la présence d'un paramètre de position, aussi, on peut utiliser cette propriété pour éviter le risque d'annulation catastrophique :
ce qui donne :
Plus est proche de la moyenne de l'échantillon, plus le calcul de la variance sera précis, mais choisir une valeur parmi les données de l'échantillon suffit à rendre le calcul stable. Si les valeurs sont petites alors il n'y aura pas de problème avec les sommes des carrés, et inversement, si elles sont grandes, alors la variance est également grandes. Dans tous les cas, le deuxième terme de la formuel est toujorus inférieur au premier, et il n'y aura donc pas d'annulation[2].
En prenant la première donnée comme paramètre , l'algorithme peut se réécrire en :
- Soient n ← 0, Sum ← 0, SumSq ← 0, K ← x0
- Pour chaque donnée x :
- n ← n + 1
- Sum ← Sum + (x – K)
- SumSq ← SumSq + (x – K) × (x – K)
- Var = (SumSq − (Sum × Sum) / n) / (n − 1)
Algorithme en deux passes
[modifier | modifier le code]Une alternative, se basant sur une autre formulation de la variance, consiste à calculer la moyenne dans un premier temps,
puis calculer la somme des carrés des différences par rapport à la moyenne,
On a donc le code :
- Soient n ← 0, Sum ← 0, SumSq ← 0
- Pour chaque donnée x :
- n ← n + 1
- Sum ← Sum + x
- Moy = Sum / n
- Pour chaque donnée x :
- SumSq ← SumSq + (x − Moy) × (x − Moy)
- Var = SumSq / (n − 1)
Cet algorithme est numériquement stable pour n petit[1],[4]. Cependant, les résultats des deux algorithmes simples (naïf et à deux passes) peuvent dépendre énormément de l'ordre des données et donc induire de mauvais résultats pour de grands ensembles de données par des erreurs répétées d'arrondi dans l'accumulation des sommes. Des techniques comme la sommation compensée peuvent aider à limiter l'impact de cette erreur.
Algorithme en ligne de Welford
[modifier | modifier le code]Il est souvent utile de calculer la variance en une seule passe et d'utiliser chaque valeur une seule fois ; par exemple, quand les données sont collectées sans espace mémoire suffisant pour stocker toutes les valeurs, ou quand le coût mémoire domine le coût calcul. Pour un tel algorithme en ligne, une relation de récurrence est nécessaire entre les données, de façon à calculer les statistiques nécessaires par une méthode stable.
Les formules suivantes peuvent être utilisées pour mettre à jour la moyenne et la variance (estimée) de la suite, pour un élément additionnel xn. Ici, désigne la moyenne de l'échantillon des n premières données , la variance de l'échantillon par un estimateur biaisé, et la variance de l'échantillon par un estimateur sans biais.
Ces formules sont sources d'instabilité numérique, car elles soustraient de manière répétée un petit nombre d'un grand nombre qui grandit avec n. Une meilleure quantité pour mettre à jour la somme des carrés des différences à la moyenne , notée ici , sera :
Cet algorithme a été trouvé par Welford[5],[6] et largement analysé[2],[7]. Il est courant de noter et [8].
- Soient n ← 0, Moy_old ← 0, Moy_new ← 0, M ← 0, S ← 0
- Pour chaque donné x :
- n ← n + 1
- M_new ← M_new + x
- Si n > 1
- M_old ← M_new
- M_new ← M_old + (x – M_old) / n
- M ← M + (x - Moy_old) × (x - Moy_new)
- Var = M / (n − 1)
Cet algorithme est moins susceptible de causer des pertes de précision par une annulation catastrophique, mais peut ne pas être aussi efficace à cause des divisions dans la boucle. Pour un algorithme en deux passes plus robuste, on peut d'abord calculer la moyenne et soustraire cette estimation des données, puis utiliser l'algorithme sur les résidus.
L'algorithme parallèle illustré infra illustre comment fusionner des ensembles multiples de statistiques calculées en ligne.
Algorithme incrémental pondéré
[modifier | modifier le code]L'algorithme peut être adapté pour prendre en compte la pondération des échantillons, en remplaçant le compteur n par la somme des poids cumulés. West (1979)[9] suggère l'algorithme suivant :
- Soient n ← 0, Sum ← 0, Sum2 ← 0, SumSq ← 0, Moy ← 0, Moy_old ← 0
- Pour chaque donnée (x , w) :
- n ← n + 1
- Sum ← Sum + w
- Moy_old ← Moy
- Moy ← Moy_old + (w / Sum) × (x – Moy_old)
- SumSq = Var + w × (x – Moy_old) × (x – Moy)
- Var = SumSq / Sum2
Algorithme parallèle
[modifier | modifier le code]Chan et al.[10] remarquent que l'algorithme en ligne de Welford décrit supra est un cas spécial d'un algorithme de somme sur deux ensembles et :
- .
Il peut être utile quand, par exemple, plusieurs unités de calcul peuvent être assignées à des parties discrètes des données d'entrée.
La méthode de Chan d'estimation de la moyenne est numériquement instable quand et sont toutes les deux très grandes, car l'erreur numérique en n'est pas échelonnée comme dan le cas où case. Dans de tels cas, on préférera .
Il peut être généralisée pour permettre la parallélisation avec AVX, avec le calcul sur carte graphique et sur clusters, et au calcul de la covariance[3].
Exemple
[modifier | modifier le code]On suppose que toutes les opérations utilisent l'arithmétique standard IEEE 754 double-precision. Pour l'échantillon (4, 7, 13, 16) d'une population infinie, la moyenne estimée est de 10, la variance estimée (sans biais) est de 30. Ces deux valeurs sont bien les résultats obtenus par l'algorithme naïf et l'algorithme à deux passes.
On considère maintenant l'échantillon (108 + 4, 108 + 7, 108 + 13, 108 + 16), qui possède la même variance que le premier. Cette fois, si l'algorithme à deux passes donne le bon résultat, l'algorithme naïf renvoie 29,333333333333332 au lieu de 30.
Tant que la perte de précision est acceptable et vue comme un défaut mineur de l'algorithme naïf, augmenter encore le décalage rendra l'erreur catastrophique. CPour l'échantillon (109 + 4, 109 + 7, 109 + 13, 109 + 16), qui est encore de variance égale à 30, l'algorithme à deux passes reste fonctionnel mais l'algorithme naïf donne −170.66666666666666. Ce problème est causé par une annulation catastrophique dans l'algorithme naïf au moment de la soustraction de deux nombres similaires à la dernière étape de l'algorithme.
Statistiques d'ordres plus élevés
[modifier | modifier le code]Terriberry[11] étend les formules de Chan pour calculer les moments centrés d'ordre 3 et 4, requis pour l'estimation de l'asymétrie et la kurtosis :
Ici les sont encore les sommes des puissances de différences par rapport à la moyenne , ce qui donne
Pour le cas incrémental (i.e., ), cela se simplifie en :
En préservant la valeur , une seule opération de division est nécessaire et les statistiques d'ordre supérieur peuvent être calculés pour des coûts incrémentaux faibles.
- Soient n ← 0, Sum ← 0, Sum2 ← 0, Sum3 ← 0, Sum4 ← 0
- Pour chaque donnée x :
- n ← n + 1
- Delta ← x – Sum
- Delta_n ← Delta / n
- Delta_2 ← Delta_n × Delta_n
- M ← Delta × Delta_n × (n – 1)
- Sum ← Sum + Delta_n
- Sum4 ← Sum4 + M × Delta_2 × (n × n – 3 × n + 3) + 6 × Delta_2 × Sum2 – 4 × Delta_n × Sum3
- Sum3 ← Sum + M × Delta_n × (n – 2) – 6 × Delta_n × Sum2
- Sum2 ← Sum2 + M
- Kurtosis ← (n × Sum4) / (Sum2 × Sum2) – 3
Pébaÿ[12] étend plus tard ces résultats aux moments centrés de tout ordre, pour les cas incrémentaux et appariés, et par la suite Pébaÿ et al.[13] pour les moments pondérés et composés. On peut également les étendre au calcul de la covariance.
Choi et Sweetman[14] offrent deux méthodes alternatives au calcul de l'asymétrie et de la kurtosis, chacune proposant des économies en coût mémoire et temps CPU pour certaines applications. La première approche est de calculer les moments statistiques en séparant les données en catégories puis en calculant les moments à partir de la géométrie de l'histogramme déduit, ce qui donne dans les faits un algorithme à une passe pour des moments d'ordre élevé. Un intérêt est que les calculs des moments statistiques peuvent être faits avec une précision arbitraire telle que les calculs peuvent être menés à la précision, par exemple, du format des stockage de données ou du matériel de mesure originel. Un histogramme relatif d'une variable aléatoire peut être construit de manière conventionnelle : l'étendue des valeurs potentielles est divisée en catégories et le nombre d'occurrences dans chaque catégorie est compté et tracé de sorte que l'aire de chaque rectangle vaut la part des valeurs d'échantillon dans chacune d'entre elles :
où et représentent la fréquence et la fréquence relative de la catégorie et est l'aire totale de l'histogramme. Après cette normalisation, les premiers moments et moments centrés de peuvent être calculés à partir de l'histogramme relatif :
où l'exposant indique les moments sont calculés à partir de l'histogramme. Pour une largeur de catégorie constante , ces deux expressions peuvent être simplifiés en utilisant :
La deuxième approche de Choi et Sweetman[14] est une méthode analytique pour combiner des moments statistiques à partir des segments individuels d'un historique tel que les moments résultants sont ceux de l'historique complet. Cette méthode peut être utilisée pour le calcul parallèle de moments statistiques avec des combinaisons déduités de ces moments, ou pour la combinaison de moments statistiques calculés de manière séquentielle.
Si ensembles de moments statistiques sont connus : pour , alors tout peut être exprimé en termes de moments équivalents :
où est généralement prise comme la durée du qe historique, ou le nombre de points si est constante.
Le bénéfice d'exprimer les moments statistiques en termes de est que les ensembles peuvent être combinés par addition, et il n'y a aucune borne supérieure sur la valeur de .
où l'indice représente l'historique concaténé ou combiné . Ces valeurs combinées de peuvent être ensuite est inversement transformées en moments bruts représentant l'historique concaténé complet
Les relations connus entre les moments bruts () et les moments centrés () sont alors utilisés pour calculer les moments centrés de l'historique concaténé. Finalement, les moments statistiques de l'historique concaténé sont calculés à partir des moments centrés :
Calcul de la covariance
[modifier | modifier le code]Des algorithmes très similaires peuvent être tirés des précédents pour le calcul de la covariance.
Algorithme naïf
[modifier | modifier le code]Un algorithme naïf de calcul de la covariance se ferait à partir de la formule :
ce qui se traduit par
- Soient n ← 0, SumX ← 0, SumY ← 0, SumXY ← 0
- Pour chaque donnée (x , y) :
- n ← n + 1
- SumX ← Sum + x
- SumY ← SumY + x
- SumXY ← SumXY + x × y
- Cov = (SumXY − (SumX × SumY) / n) / (n − 1)
Avec estimation de la moyenne
[modifier | modifier le code]Comme pour la variance, la covariance de deux vecteurs aléatoires est invariante par translation, donc pour toutes valeurs constantes et données, on a :
et choisir une valeur parmi les données va rendre le calcul plus stable contre les annulations catastrophiques et plus robuste contre les grandes sommes.
À deux passes
[modifier | modifier le code]L'algorithme à deux passes calcul d'abord les deux moyennes, puis la covariance :
On peut écrire l'algorithme comme suit :
- Soient n ← 0, SumX ← 0, SumY ← 0, SumXY ← 0
- Pour chaque donnée (x ,y) :
- n ← n + 1
- SumX ← SumX + x
- SumY ← SumY + y
- MoyX ← SumX / n
- MoyY ← SumY / n
- Pour chaque donnée (x ,y) :
- SumXY ← SumXY + (x – MoyX) × (y – MoyY)
- Cov = SumXY / n
Une version compensée légèrement plus précise exécute l'algorithme naïf complet sur les résidus. Les sommes finales et doivent être nulles, mais la deuxième passe compense toute petite erreur.
En ligne
[modifier | modifier le code]Un algorithme à une passe stable, similaire à l'algorithme en ligne pour le calcul de la variance, qui calcule le co-moment :
L'asymétrie apparente dans la dernière équation est due au fait que , donc les deux termes ajoutés sont égaux à . Même une précision plus grande peut être atteinte en calculant d'abord les moyennes, puis en utilisant l'algorithme stable à une passe sur les résidus.
Alors la covariance peut être calculé par
De même, il y a une formule pour combiner les covariances de deux ensembles qui peut être utilisé pour paralléliser le calcul[3]:
Version groupée pondérée
[modifier | modifier le code]Une version de l'algorithme en ligne pondérée qui met à jour par groupe existe aussi : soient les poids, et on écrit
La covariance s'obtient avec
Voir aussi
[modifier | modifier le code]Références
[modifier | modifier le code]- (en) Cet article est partiellement ou en totalité issu de l’article de Wikipédia en anglais intitulé « Algorithms for calculating variance » (voir la liste des auteurs).
- Bo Einarsson, Accuracy and Reliability in Scientific Computing, SIAM, (ISBN 978-0-89871-584-2, lire en ligne), p. 47
- (en) Tony F. Chan, Gene H. Golub et Randall J. LeVeque, « Algorithms for computing the sample variance: Analysis and recommendations », The American Statistician, vol. 37, no 3, , p. 242–247 (DOI 10.1080/00031305.1983.10483115, JSTOR 2683386, lire en ligne [archive du ])
- (en) Erich Schubert et Michael Gertz, Numerically stable parallel computation of (co-)variance, ACM, , 10 p. (ISBN 9781450365055, DOI 10.1145/3221269.3223036, S2CID 49665540, lire en ligne)
- (en) Nicholas J. Higham, Accuracy and Stability of Numerical Algorithms, Philadelphia, PA, Society for Industrial and Applied Mathematics, (ISBN 978-0-89871-802-7, DOI 10.1137/1.9780898718027, lire en ligne), « Problem 1.10 ». Metadata le liste dans ACM Digital Library.
- (en) B. P. Welford, « Note on a method for calculating corrected sums of squares and products », Technometrics, vol. 4, no 3, , p. 419–420 (DOI 10.2307/1266577, JSTOR 1266577)
- Donald E. Knuth (1998). The Art of Computer Programming, volume 2: Seminumerical Algorithms, 3rd edn., p. 232. Boston: Addison-Wesley.
- (en) Robert F. Ling, « Comparison of Several Algorithms for Computing Sample Means and Variances », Journal of the American Statistical Association, vol. 69, no 348, , p. 859–866 (DOI 10.2307/2286154, JSTOR 2286154)
- (en) John D. Cook, « Accurately computing sample variance », sur John D. Cook Consulting: Expert consulting in applied mathematics & data privacy,
- (en) D. H. D. West, « Updating Mean and Variance Estimates: An Improved Method », Communications of the ACM, vol. 22, no 9, , p. 532–535 (DOI 10.1145/359146.359153
, S2CID 30671293)
- (en) Tony F. Chan, Gene H. Golub et Randall J. LeVeque, « Updating Formulae and a Pairwise Algorithm for Computing Sample Variances », Department of Computer Science, Stanford University,
- (en) Timothy B. Terriberry, « Computing Higher-Order Moments Online » [archive du ], (consulté le )
- (en) Philippe Pierre Pébay (Sponsoring Org.: USDOE), « Formulas for Robust, One-Pass Parallel Computation of Covariances and Arbitrary-Order Statistical Moments », sur UNT Digital Library, Albuquerque, NM, and Livermore, CA (United States), Sandia National Laboratories (SNL), (DOI 10.2172/1028931)
- (en) Philippe Pébaÿ, Timothy Terriberry, Hemanth Kolla et Janine Bennett, « Numerically Stable, Scalable Formulas for Parallel and Online Computation of Higher-Order Multivariate Central Moments with Arbitrary Weights », Springer, vol. 31, no 4, , p. 1305–1325 (DOI 10.1007/s00180-015-0637-z, S2CID 124570169, lire en ligne)
- (en) Myoungkeun Choi et Bert Sweetman, « Efficient Calculation of Statistical Moments for Structural Health Monitoring », Journal of Structural Health Monitoring, vol. 9, no 1, , p. 13–24 (DOI 10.1177/1475921709341014, S2CID 17534100)
Liens externes
[modifier | modifier le code]- (en) Eric W. Weisstein, « Sample Variance Computation », sur MathWorld