# 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

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, la formule est :
    $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, on peut proposer un algorithme naïf pour calculer une estimation de la variance :
* 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$^,$, aussi plusieurs algorithmes alternatifs, numériquement stables, ont été proposés. 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 grande. Dans tous les cas, le deuxième terme de la formule est toujours inférieur au premier, et il n'y aura donc pas d'annulation.
En prenant la première donnée comme paramètre $K$, l'algorithme peut se réécrire en :
* Soient n ← 0, Sum ← 0, SumSq ← 0, K ← x$_0$
* 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

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 ← 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$^,$. 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, x$^¯$$_n$=1/n ∑$_{i=1}$$^n$x$_i${\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})$, σ$_n$$^2$=1/n ∑$_{i=1}$$^n$(x$_i$−x$^¯$$_n$)$^2${\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 s$_n$$^2$=1/n−1 ∑$_{i=1}$$^n$(x$_i$−x$^¯$$_n$)$^2${\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 ∑$_{i=1}$$^n$(x$_i$−x$^¯$$_n$)$^2${\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$^,$ et largement analysé$^,$. Il est courant de noter $M_{k}={\bar {x}}_{k}$ et $S_{k}=M_{2,k}$.
* 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é

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

Chan et al. 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 tous les deux très grands, car l'erreur numérique en $\delta ={\bar {x}}_{B}-{\bar {x}}_{A}$ n'est pas échelonnée comme dans le cas où $n_{B}=1$. Dans de tels cas, on préférera x$^¯$$_{AB}$=(n$_A$x$^¯$$_A$+n$_B$x$^¯$$_B$)/n$_{AB}$ {\textstyle {\bar {x}}_{AB}={\frac {n_{A}{\bar {x}}_{A}+n_{B}{\bar {x}}_{B}}{n_{AB}}}}.
Cela peut être généralisé pour permettre la parallélisation avec AVX, avec le calcul sur carte graphique et sur clusters, et au calcul de la covariance.

## 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$^8$ + 4, 10$^8$ + 7, 10$^8$ + 13, 10$^8$ + 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. Pour l'échantillon (10$^9$ + 4, 10$^9$ + 7, 10$^9$ + 13, 10$^9$ + 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 é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 ∑(x−x$^¯$)$^k${\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ées 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ÿ é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. pour les moments pondérés et composés. On peut également les étendre au calcul de la covariance.
Choi et Sweetman 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 de 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 A=∑$_{k=1}$$^K$h(x$_k$)Δx$_k${\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 calculés à partir de l'histogramme. Pour une largeur de catégorie constante $\Delta x_{k}=\Delta x$, ces deux expressions peuvent être simplifiées 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 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éduites 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 pris comme la durée du qe 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 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 connues entre les moments bruts ($m_{n}$) et les moments centrés ($\theta _{n}=\operatorname {E} [(x-\mu )^{n}])$) sont alors utilisées 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 ← 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

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 calcule 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 ← 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 ∑$_i$x$_i${\textstyle \sum _{i}x_{i}} et ∑$_i$y$_i${\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 C$_n$=∑$_{i=1}$$^n$(x$_i$−x$^¯$$_n$)(y$_i$−y$^¯$$_n$){\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 (x$_n$−x$^¯$$_n$)=n−1/n (x$_n$−x$^¯$$_{n−1}$){\textstyle (x_{n}-{\bar {x}}_{n})={\frac {n-1}{n}}(x_{n}-{\bar {x}}_{n-1})}, donc les deux termes ajoutés sont égaux à n−1/n (x$_n$−x$^¯$$_{n−1}$)(y$_n$−y$^¯$$_{n−1}$){\textstyle {\frac {n-1}{n}}(x_{n}-{\bar {x}}_{n-1})(y_{n}-{\bar {y}}_{n-1})}. Une précision plus grande peut même ê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ée 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ée pour paralléliser le calcul:
    $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}}}$