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

Si ce bandeau n'est plus pertinent, retirez-le. Cliquez ici pour en savoir plus.

Cette section ou cet article est une traduction incomplète (juin 2024).

Vous pouvez modifier la page pour effectuer la traduction.

Algorithme naïf

Une formule usuelle pour le calcul de la variance d'une population de taille N est :

\sigma ^{2}={\overline {(x^{2})}}-{\bar {x}}^{2}={\frac {1}{N}}\left[\sum _{i=1}^{N}x_{i}^{2}-{\frac {1}{N}}\left(\sum _{i=1}^{N}x_{i}\right)^{2}\right].

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:

s^{2}={\frac {n}{n-1}}\left({\frac {\sum _{i=1}^{n}x_{i}^{2}}{n}}-\left({\frac {\sum _{i=1}^{n}x_{i}}{n}}\right)^{2}\right).

Ainsi, un algorithme naïve algorithm to calculate the estimated variance is given by the following:

Soient $n \leftarrow 0, Sum \leftarrow 0, SumSq \leftarrow 0$
Pour chaque donné x :
- $n \leftarrow n + 1$
- $Sum \leftarrow Sum + x$
- $SumSq \leftarrow SumSq + x \times x$
$Var = (SumSq - (Sum \times 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 \times 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

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 :

\forall K\in \mathbb {R} ,\operatorname {Var} (X-K)=\operatorname {Var} (X).

ce qui donne :

\sigma ^{2}={\frac {\sum _{i=1}^{n}x_{i}^{2}-(\sum _{i=1}^{n}x_{i})^{2}/n}{n-1}}={\frac {\sum _{i=1}^{n}(x_{i}-K)^{2}-(\sum _{i=1}^{n}(x_{i}-K))^{2}/n}{n-1}}.

Plus $K$ 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 $(x_{i}-K)$ 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 $K$ , l'algorithme peut se réécrire en :

Soient $n \leftarrow 0, Sum \leftarrow 0, SumSq \leftarrow 0, K \leftarrow x 0$
Pour chaque donnée x :
- $n \leftarrow n + 1$
- $Sum \leftarrow Sum + (x - K)$
- $SumSq \leftarrow SumSq + (x - K) \times (x - K)$
$Var = (SumSq - (Sum \times Sum) / n) / (n - 1)$

Algorithme en deux passes

Une alternative, se basant sur une autre formulation de la variance, consiste à calculer la moyenne dans un premier temps,

{\bar {x}}={\frac {\sum _{j=1}^{n}x_{j}}{n}},

puis calculer la somme des carrés des différences par rapport à la moyenne,

s^{2}={\dfrac {\sum _{i=1}^{n}(x_{i}-{\bar {x}})^{2}}{n-1}}.

On a donc le code :

Soient $n \leftarrow 0, Sum \leftarrow 0, SumSq \leftarrow 0$
Pour chaque donnée x :
- $n \leftarrow n + 1$
- $Sum \leftarrow Sum + x$
$Moy = Sum / n$
Pour chaque donnée x :
- $SumSq \leftarrow SumSq + (x - Moy) \times (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

Il est souvent utile de calculer la variance en une seule passe et d'utiliser chaque valeur $x_{i}$ 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 x_n. Ici, ${\textstyle {\overline {x}}_{n}={\frac {1}{n}}\sum _{i=1}^{n}x_{i}}$ désigne la moyenne de l'échantillon des n premières données $(x_{1},\dots ,x_{n})$ , ${\textstyle \sigma _{n}^{2}={\frac {1}{n}}\sum _{i=1}^{n}\left(x_{i}-{\overline {x}}_{n}\right)^{2}}$ la variance de l'échantillon par un estimateur biaisé, et ${\textstyle s_{n}^{2}={\frac {1}{n-1}}\sum _{i=1}^{n}\left(x_{i}-{\overline {x}}_{n}\right)^{2}}$ la variance de l'échantillon par un estimateur sans biais.

{\bar {x}}_{n}={\frac {(n-1)\,{\bar {x}}_{n-1}+x_{n}}{n}}={\bar {x}}_{n-1}+{\frac {x_{n}-{\bar {x}}_{n-1}}{n}}

\sigma _{n}^{2}={\frac {(n-1)\,\sigma _{n-1}^{2}+(x_{n}-{\bar {x}}_{n-1})(x_{n}-{\bar {x}}_{n})}{n}}=\sigma _{n-1}^{2}+{\frac {(x_{n}-{\bar {x}}_{n-1})(x_{n}-{\bar {x}}_{n})-\sigma _{n-1}^{2}}{n}}.

s_{n}^{2}={\frac {n-2}{n-1}}\,s_{n-1}^{2}+{\frac {(x_{n}-{\bar {x}}_{n-1})^{2}}{n}}=s_{n-1}^{2}+{\frac {(x_{n}-{\bar {x}}_{n-1})^{2}}{n}}-{\frac {s_{n-1}^{2}}{n-1}},\quad n>1

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 ${\textstyle \sum _{i=1}^{n}(x_{i}-{\bar {x}}_{n})^{2}}$ , notée ici $M_{2,n}$ , sera :

{\begin{aligned}M_{2,n}&=M_{2,n-1}+(x_{n}-{\bar {x}}_{n-1})(x_{n}-{\bar {x}}_{n})\\[4pt]\sigma _{n}^{2}&={\frac {M_{2,n}}{n}}\\[4pt]s_{n}^{2}&={\frac {M_{2,n}}{n-1}}\end{aligned}}

Cet algorithme a été trouvé par Welford^[5]^,^[6] et largement analysé^[2]^,^[7]. Il est courant de noter $M_{k}={\bar {x}}_{k}$ et $S_{k}=M_{2,k}$ ^[8].

Soient $n \leftarrow 0, Moy_old \leftarrow 0, Moy_new \leftarrow 0, M \leftarrow 0, S \leftarrow 0$
Pour chaque donné x :
- $n \leftarrow n + 1$
- $M_new \leftarrow M_new + x$
- Si n > 1
  - $M_old \leftarrow M_new$
  - $M_new \leftarrow M_old + (x - M_old) / n$
  - $M \leftarrow M + (x - Moy_old) \times (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é

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 \leftarrow 0, Sum \leftarrow 0, Sum2 \leftarrow 0, SumSq \leftarrow 0, Moy \leftarrow 0, Moy_old \leftarrow 0$
Pour chaque donnée (x , w) :
- $n \leftarrow n + 1$
- $Sum \leftarrow Sum + w$
- $Moy_old \leftarrow Moy$
- $Moy \leftarrow Moy_old + (w / Sum) \times (x - Moy_old)$
- $SumSq = Var + w \times (x - Moy_old) \times (x - Moy)$
$Var = SumSq / Sum2$

Algorithme parallèle

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 $A$ et $B$ :

{\begin{aligned}n_{AB}&=n_{A}+n_{B}\\\delta &={\bar {x}}_{B}-{\bar {x}}_{A}\\{\bar {x}}_{AB}&={\bar {x}}_{A}+\delta \cdot {\frac {n_{B}}{n_{AB}}}\\M_{2,AB}&=M_{2,A}+M_{2,B}+\delta ^{2}\cdot {\frac {n_{A}n_{B}}{n_{AB}}}\\\end{aligned}}

.

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 $n_{A}\approx n_{B}$ et sont toutes les deux très grandes, car l'erreur numérique en $\delta ={\bar {x}}_{B}-{\bar {x}}_{A}$ n'est pas échelonnée comme dan le cas où $n_{B}=1$ case. Dans de tels cas, on préférera ${\textstyle {\bar {x}}_{AB}={\frac {n_{A}{\bar {x}}_{A}+n_{B}{\bar {x}}_{B}}{n_{AB}}}}$ .

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

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 (10⁸ + 4, 10⁸ + 7, 10⁸ + 13, 10⁸ + 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 (10⁹ + 4, 10⁹ + 7, 10⁹ + 13, 10⁹ + 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

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 :

{\begin{aligned}M_{3,X}=M_{3,A}+M_{3,B}&{}+\delta ^{3}{\frac {n_{A}n_{B}(n_{A}-n_{B})}{n_{X}^{2}}}+3\delta {\frac {n_{A}M_{2,B}-n_{B}M_{2,A}}{n_{X}}}\\[6pt]M_{4,X}=M_{4,A}+M_{4,B}&{}+\delta ^{4}{\frac {n_{A}n_{B}\left(n_{A}^{2}-n_{A}n_{B}+n_{B}^{2}\right)}{n_{X}^{3}}}\\[6pt]&{}+6\delta ^{2}{\frac {n_{A}^{2}M_{2,B}+n_{B}^{2}M_{2,A}}{n_{X}^{2}}}+4\delta {\frac {n_{A}M_{3,B}-n_{B}M_{3,A}}{n_{X}}}\end{aligned}}

Ici les $M_{k}$ sont encore les sommes des puissances de différences par rapport à la moyenne ${\textstyle \sum (x-{\overline {x}})^{k}}$ , ce qui donne

{\begin{aligned}&{\text{skewness}}=g_{1}={\frac {{\sqrt {n}}M_{3}}{M_{2}^{3/2}}},\\[4pt]&{\text{kurtosis}}=g_{2}={\frac {nM_{4}}{M_{2}^{2}}}-3.\end{aligned}}

Pour le cas incrémental (i.e., $B=\{x\}$ ), cela se simplifie en :

{\begin{aligned}\delta &=x-m\\[5pt]m'&=m+{\frac {\delta }{n}}\\[5pt]M_{2}'&=M_{2}+\delta ^{2}{\frac {n-1}{n}}\\[5pt]M_{3}'&=M_{3}+\delta ^{3}{\frac {(n-1)(n-2)}{n^{2}}}-{\frac {3\delta M_{2}}{n}}\\[5pt]M_{4}'&=M_{4}+{\frac {\delta ^{4}(n-1)(n^{2}-3n+3)}{n^{3}}}+{\frac {6\delta ^{2}M_{2}}{n^{2}}}-{\frac {4\delta M_{3}}{n}}\end{aligned}}

En préservant la valeur $\delta /n$ , 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 \leftarrow 0, Sum \leftarrow 0, Sum2 \leftarrow 0, Sum3 \leftarrow 0, Sum4 \leftarrow 0$
Pour chaque donnée x :
- $n \leftarrow n + 1$
- $Delta \leftarrow x - Sum$
- $Delta_n \leftarrow Delta / n$
- $Delta_2 \leftarrow Delta_n \times Delta_n$
- $M \leftarrow Delta \times Delta_n \times (n - 1)$
- $Sum \leftarrow Sum + Delta_n$
- $Sum4 \leftarrow Sum4 + M \times Delta_2 \times (n \times n - 3 \times n + 3) + 6 \times Delta_2 \times Sum2 - 4 \times Delta_n \times Sum3$
- $Sum3 \leftarrow Sum + M \times Delta_n \times (n - 2) - 6 \times Delta_n \times Sum2$
- $Sum2 \leftarrow Sum2 + M$
$Kurtosis \leftarrow (n \times Sum4) / (Sum2 \times 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 :

H(x_{k})={\frac {h(x_{k})}{A}}

où $h(x_{k})$ et $H(x_{k})$ représentent la fréquence et la fréquence relative de la catégorie $x_{k}$ et ${\textstyle A=\sum _{k=1}^{K}h(x_{k})\,\Delta x_{k}}$ est l'aire totale de l'histogramme. Après cette normalisation, les $n$ premiers moments et moments centrés de $x(t)$ peuvent être calculés à partir de l'histogramme relatif :

m_{n}^{(h)}=\sum _{k=1}^{K}x_{k}^{n}H(x_{k})\,\Delta x_{k}={\frac {1}{A}}\sum _{k=1}^{K}x_{k}^{n}h(x_{k})\,\Delta x_{k}

\theta _{n}^{(h)}=\sum _{k=1}^{K}\left(x_{k}-m_{1}^{(h)}\right)^{n}\,H(x_{k})\,\Delta x_{k}={\frac {1}{A}}\sum _{k=1}^{K}\left(x_{k}-m_{1}^{(h)}\right)^{n}h(x_{k})\,\Delta x_{k}

où l'exposant $^{(h)}$ indique les moments sont calculés à partir de l'histogramme. Pour une largeur de catégorie constante $\Delta x_{k}=\Delta x$ , ces deux expressions peuvent être simplifiés en utilisant $I=A/\Delta x$ :

m_{n}^{(h)}={\frac {1}{I}}\sum _{k=1}^{K}x_{k}^{n}\,h(x_{k})

\theta _{n}^{(h)}={\frac {1}{I}}\sum _{k=1}^{K}\left(x_{k}-m_{1}^{(h)}\right)^{n}h(x_{k})

La deuxième approche de Choi et Sweetman^[14] est une méthodologie 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éthodologie peut être utilisé 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 $Q$ ensembles de moments statistiques sont connus : $(\gamma _{0,q},\mu _{q},\sigma _{q}^{2},\alpha _{3,q},\alpha _{4,q})\quad$ pour $q=1,2,\ldots ,Q$ , alors tout $\gamma _{n}$ peut être exprimé en termes de $n$ moments équivalents :

\gamma _{n,q}=m_{n,q}\gamma _{0,q}\qquad \quad {\textrm {pour}}\quad n=1,2,3,4\quad {\text{ et }}\quad q=1,2,\dots ,Q

où $\gamma _{0,q}$ est généralement prise comme la durée du $q e$ historique, ou le nombre de points si $\Delta t$ est constante.

Le bénéfice d'exprimer les moments statistiques en termes de $\gamma$ est que les $Q$ ensembles peuvent être combinés par addition, et il n'y a aucune borne supérieure sur la valeur de $Q$ .

\gamma _{n,c}=\sum _{q=1}^{Q}\gamma _{n,q}\quad \quad {\text{pour }}n=0,1,2,3,4

où l'indice $_{c}$ représente l'historique concaténé ou combiné $\gamma$ . Ces valeurs combinées de $\gamma$ peuvent être ensuite est inversement transformées en moments bruts représentant l'historique concaténé complet

m_{n,c}={\frac {\gamma _{n,c}}{\gamma _{0,c}}}\quad {\text{for }}n=1,2,3,4

Les relations connus entre les moments bruts ( $m_{n}$ ) et les moments centrés ( $\theta _{n}=\operatorname {E} [(x-\mu )^{n}])$ ) 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 :

\mu _{c}=m_{1,c}\qquad \sigma _{c}^{2}=\theta _{2,c}\qquad \alpha _{3,c}={\frac {\theta _{3,c}}{\sigma _{c}^{3}}}\qquad \alpha _{4,c}={\frac {\theta _{4,c}}{\sigma _{c}^{4}}}-3

Calcul de la covariance

Des algorithmes très similaires peuvent être tirés des précédents pour le calcul de la covariance.

Algorithme naïf

Un algorithme naïf de calcul de la covariance se ferait à partir de la formule :

\operatorname {Cov} (X,Y)={\frac {\sum _{i=1}^{n}x_{i}y_{i}-(\sum _{i=1}^{n}x_{i})(\sum _{i=1}^{n}y_{i})/n}{n}}.

ce qui se traduit par

Soient $n \leftarrow 0, SumX \leftarrow 0, SumY \leftarrow 0, SumXY \leftarrow 0$
Pour chaque donnée (x , y) :
- $n \leftarrow n + 1$
- $SumX \leftarrow Sum + x$
- $SumY \leftarrow SumY + x$
- $SumXY \leftarrow SumXY + x \times y$
$Cov = (SumXY - (SumX \times SumY) / n) / (n - 1)$

Avec estimation de la moyenne

Comme pour la variance, la covariance de deux vecteurs aléatoires est invariante par translation, donc pour toutes valeurs constantes $k_{x}$ et $k_{y}$ données, on a :

\operatorname {Cov} (X,Y)=\operatorname {Cov} (X-k_{x},Y-k_{y})={\dfrac {\sum _{i=1}^{n}(x_{i}-k_{x})(y_{i}-k_{y})-(\sum _{i=1}^{n}(x_{i}-k_{x}))(\sum _{i=1}^{n}(y_{i}-k_{y}))/n}{n}}.

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

L'algorithme à deux passes calcul d'abord les deux moyennes, puis la covariance :

{\bar {x}}=\sum _{i=1}^{n}x_{i}/n\ ,\ {\bar {y}}=\sum _{i=1}^{n}y_{i}/n

\operatorname {Cov} (X,Y)={\frac {\sum _{i=1}^{n}(x_{i}-{\bar {x}})(y_{i}-{\bar {y}})}{n}}.

On peut écrire l'algorithme comme suit :

Soient $n \leftarrow 0, SumX \leftarrow 0, SumY \leftarrow 0, SumXY \leftarrow 0$
Pour chaque donnée (x ,y) :
- $n \leftarrow n + 1$
- $SumX \leftarrow SumX + x$
- $SumY \leftarrow SumY + y$
$MoyX \leftarrow SumX / n$
$MoyY \leftarrow SumY / n$
Pour chaque donnée (x ,y) :
- $SumXY \leftarrow SumXY + (x - MoyX) \times (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 ${\textstyle \sum _{i}x_{i}}$ et ${\textstyle \sum _{i}y_{i}}$ doivent être nulles, mais la deuxième passe compense toute petite erreur.

En ligne

Un algorithme à une passe stable, similaire à l'algorithme en ligne pour le calcul de la variance, qui calcule le co-moment ${\textstyle C_{n}=\sum _{i=1}^{n}(x_{i}-{\bar {x}}_{n})(y_{i}-{\bar {y}}_{n})}$ :

{\begin{alignedat}{2}{\bar {x}}_{n}&={\bar {x}}_{n-1}&\,+\,&{\frac {x_{n}-{\bar {x}}_{n-1}}{n}}\\[5pt]{\bar {y}}_{n}&={\bar {y}}_{n-1}&\,+\,&{\frac {y_{n}-{\bar {y}}_{n-1}}{n}}\\[5pt]C_{n}&=C_{n-1}&\,+\,&(x_{n}-{\bar {x}}_{n})(y_{n}-{\bar {y}}_{n-1})\\[5pt]&=C_{n-1}&\,+\,&(x_{n}-{\bar {x}}_{n-1})(y_{n}-{\bar {y}}_{n})\end{alignedat}}

L'asymétrie apparente dans la dernière équation est due au fait que ${\textstyle (x_{n}-{\bar {x}}_{n})={\frac {n-1}{n}}(x_{n}-{\bar {x}}_{n-1})}$ , donc les deux termes ajoutés sont égaux à ${\textstyle {\frac {n-1}{n}}(x_{n}-{\bar {x}}_{n-1})(y_{n}-{\bar {y}}_{n-1})}$ . 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

{\begin{aligned}\operatorname {Cov} _{N}(X,Y)={\frac {C_{N}}{N}}&={\frac {\operatorname {Cov} _{N-1}(X,Y)\cdot (N-1)+(x_{n}-{\bar {x}}_{n})(y_{n}-{\bar {y}}_{n-1})}{N}}\\&={\frac {\operatorname {Cov} _{N-1}(X,Y)\cdot (N-1)+(x_{n}-{\bar {x}}_{n-1})(y_{n}-{\bar {y}}_{n})}{N}}\\&={\frac {\operatorname {Cov} _{N-1}(X,Y)\cdot (N-1)+{\frac {N-1}{N}}(x_{n}-{\bar {x}}_{n-1})(y_{n}-{\bar {y}}_{n-1})}{N}}\\&={\frac {\operatorname {Cov} _{N-1}(X,Y)\cdot (N-1)+{\frac {N}{N-1}}(x_{n}-{\bar {x}}_{n})(y_{n}-{\bar {y}}_{n})}{N}}.\end{aligned}}

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

C_{X}=C_{A}+C_{B}+({\bar {x}}_{A}-{\bar {x}}_{B})({\bar {y}}_{A}-{\bar {y}}_{B})\cdot {\frac {n_{A}n_{B}}{n_{X}}}.

Version groupée pondérée

Une version de l'algorithme en ligne pondérée qui met à jour par groupe existe aussi : soient $w_{1},\dots w_{N}$ les poids, et on écrit

{\begin{alignedat}{2}{\bar {x}}_{n+k}&={\bar {x}}_{n}&\,+\,&{\frac {\sum _{i=n+1}^{n+k}w_{i}(x_{i}-{\bar {x}}_{n})}{\sum _{i=1}^{n+k}w_{i}}}\\{\bar {y}}_{n+k}&={\bar {y}}_{n}&\,+\,&{\frac {\sum _{i=n+1}^{n+k}w_{i}(y_{i}-{\bar {y}}_{n})}{\sum _{i=1}^{n+k}w_{i}}}\\C_{n+k}&=C_{n}&\,+\,&\sum _{i=n+1}^{n+k}w_{i}(x_{i}-{\bar {x}}_{n+k})(y_{i}-{\bar {y}}_{n})\\&=C_{n}&\,+\,&\sum _{i=n+1}^{n+k}w_{i}(x_{i}-{\bar {x}}_{n})(y_{i}-{\bar {y}}_{n+k})\\\end{alignedat}}

La covariance s'obtient avec

\operatorname {Cov} _{N}(X,Y)={\frac {C_{N}}{\sum _{i=1}^{N}w_{i}}}

Voir aussi

Références

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

↑ ^{a et b} Bo Einarsson, Accuracy and Reliability in Scientific Computing, SIAM, 2005 (ISBN 978-0-89871-584-2, lire en ligne), p. 47
↑ ^{a b et c} (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, n^o 3,‎ 1983, p. 242–247 (DOI 10.1080/00031305.1983.10483115, JSTOR 2683386, lire en ligne [archive du 9 octobre 2022])
↑ ^{a b et c} (en) Erich Schubert et Michael Gertz, Numerically stable parallel computation of (co-)variance, ACM, 9 juillet 2018, 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, 2002 (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, n^o 3,‎ 1962, 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, n^o 348,‎ 1974, 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, 30 septembre 2022
↑ (en) D. H. D. West, « Updating Mean and Variance Estimates: An Improved Method », Communications of the ACM, vol. 22, n^o 9,‎ 1979, 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, novembre 1979
↑ (en) Timothy B. Terriberry, « Computing Higher-Order Moments Online » [archive du 23 avril 2014], 15 octobre 2008 (consulté le 5 mai 2008)
↑ (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), septembre 2008 (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, n^o 4,‎ 2016, p. 1305–1325 (DOI 10.1007/s00180-015-0637-z, S2CID 124570169, lire en ligne)
↑ ^{a et b} (en) Myoungkeun Choi et Bert Sweetman, « Efficient Calculation of Statistical Moments for Structural Health Monitoring », Journal of Structural Health Monitoring, vol. 9, n^o 1,‎ 2010, p. 13–24 (DOI 10.1177/1475921709341014, S2CID 17534100)

Liens externes

(en) Eric W. Weisstein, « Sample Variance Computation », sur MathWorld

Portail des probabilités et de la statistique

[Einarsson2005-1] {a et b} Bo Einarsson, Accuracy and Reliability in Scientific Computing, SIAM, 2005 (ISBN 978-0-89871-584-2, lire en ligne), p. 47

[Chan1983-2] {a b et c} (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, n^o 3,‎ 1983, p. 242–247 (DOI 10.1080/00031305.1983.10483115, JSTOR 2683386, lire en ligne [archive du 9 octobre 2022])

[:1-3] {a b et c} (en) Erich Schubert et Michael Gertz, Numerically stable parallel computation of (co-)variance, ACM, 9 juillet 2018, 10 p. (ISBN 9781450365055, DOI 10.1145/3221269.3223036, S2CID 49665540, lire en ligne)

[4] (en) Nicholas J. Higham, Accuracy and Stability of Numerical Algorithms, Philadelphia, PA, Society for Industrial and Applied Mathematics, 2002 (ISBN 978-0-89871-802-7, DOI 10.1137/1.9780898718027, lire en ligne), « Problem 1.10 ». Metadata le liste dans ACM Digital Library.

[5] (en) B. P. Welford, « Note on a method for calculating corrected sums of squares and products », Technometrics, vol. 4, n^o 3,‎ 1962, p. 419–420 (DOI 10.2307/1266577, JSTOR 1266577)

[6] Donald E. Knuth (1998). The Art of Computer Programming, volume 2: Seminumerical Algorithms, 3rd edn., p. 232. Boston: Addison-Wesley.

[7] (en) Robert F. Ling, « Comparison of Several Algorithms for Computing Sample Means and Variances », Journal of the American Statistical Association, vol. 69, n^o 348,‎ 1974, p. 859–866 (DOI 10.2307/2286154, JSTOR 2286154)

[8] (en) John D. Cook, « Accurately computing sample variance », sur John D. Cook Consulting: Expert consulting in applied mathematics & data privacy, 30 septembre 2022

[9] (en) D. H. D. West, « Updating Mean and Variance Estimates: An Improved Method », Communications of the ACM, vol. 22, n^o 9,‎ 1979, p. 532–535 (DOI 10.1145/359146.359153 , S2CID 30671293)

[:0-10] (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, novembre 1979

[11] (en) Timothy B. Terriberry, « Computing Higher-Order Moments Online » [archive du 23 avril 2014], 15 octobre 2008 (consulté le 5 mai 2008)

[12] (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), septembre 2008 (DOI 10.2172/1028931)

[13] (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, n^o 4,‎ 2016, p. 1305–1325 (DOI 10.1007/s00180-015-0637-z, S2CID 124570169, lire en ligne)

[Choi2010-14] {a et b} (en) Myoungkeun Choi et Bert Sweetman, « Efficient Calculation of Statistical Moments for Structural Health Monitoring », Journal of Structural Health Monitoring, vol. 9, n^o 1,‎ 2010, p. 13–24 (DOI 10.1177/1475921709341014, S2CID 17534100)

[1]

[2]

[3]

[4]

[5]

[6]

[7]

[8]

[9]

[10]

[11]

[12]

[13]

[14]