Utilisateur:Kelam/Bac à sable perso

Ici, on se présente

Ici, on cause

Ici, on ébauche

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, 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 ← 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. This is particularly bad if the standard deviation is small relative to the mean.

Calcul avec des données décalées

La variance est invariant en présence d'un paramètre de position, aussi, on peut utiliser cette propriété pour éviter le risque d'annulation catastrophique :

$\forall Kin\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.

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₀
• 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
• n ← 0
• Pour chaque donnée x :
  • n ← n + 1
  • 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, ${\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 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 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.

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.

Higher-order statistics

Terriberry extends Chan's formulae to calculating the third and fourth central moments, needed for example when estimating skewness and 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}}$

Here the $M_{k}$ are again the sums of powers of differences from the mean ${\textstyle \sum (x-{\overline {x}})^{k}}$, giving

${\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}}$

For the incremental case (i.e., $B=\{x\}$), this simplifies to:

${\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}}$

By preserving the value $\delta /n$, only one division operation is needed and the higher-order statistics can thus be calculated for little incremental cost.

An example of the online algorithm for kurtosis implemented as described is:

def online_kurtosis(data):
    n = mean = M2 = M3 = M4 = 0

    for x in data:
        n1 = n
        n = n + 1
        delta = x - mean
        delta_n = delta / n
        delta_n2 = delta_n**2
        term1 = delta * delta_n * n1
        mean = mean + delta_n
        M4 = M4 + term1 * delta_n2 * (n**2 - 3*n + 3) + 6 * delta_n2 * M2 - 4 * delta_n * M3
        M3 = M3 + term1 * delta_n * (n - 2) - 3 * delta_n * M2
        M2 = M2 + term1

    # Note, you may also calculate variance using M2, and skewness using M3
    # Caution: If all the inputs are the same, M2 will be 0, resulting in a division by 0.
    kurtosis = (n * M4) / (M2**2) - 3
    return kurtosis

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

Pébaÿ further extends these results to arbitrary-order central moments, for the incremental and the pairwise cases, and subsequently Pébaÿ et al. for weighted and compound moments. One can also find there similar formulas for covariance.

Choi and Sweetman offer two alternative methods to compute the skewness and kurtosis, each of which can save substantial computer memory requirements and CPU time in certain applications. The first approach is to compute the statistical moments by separating the data into bins and then computing the moments from the geometry of the resulting histogram, which effectively becomes a one-pass algorithm for higher moments. One benefit is that the statistical moment calculations can be carried out to arbitrary accuracy such that the computations can be tuned to the precision of, e.g., the data storage format or the original measurement hardware. A relative histogram of a random variable can be constructed in the conventional way: the range of potential values is divided into bins and the number of occurrences within each bin are counted and plotted such that the area of each rectangle equals the portion of the sample values within that bin:

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

where $h(x_{k})$ and $H(x_{k})$ represent the frequency and the relative frequency at bin $x_{k}$ and ${\textstyle A=\sum _{k=1}^{K}h(x_{k})\,\Delta x_{k}}$ is the total area of the histogram. After this normalization, the $n$ raw moments and central moments of $x(t)$ can be calculated from the relative histogram:

$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}{\Big (}x_{k}-m_{1}^{(h)}{\Big )}^{n}\,H(x_{k})\,\Delta x_{k}={\frac {1}{A}}\sum _{k=1}^{K}{\Big (}x_{k}-m_{1}^{(h)}{\Big )}^{n}h(x_{k})\,\Delta x_{k}$

where the superscript $^{(h)}$ indicates the moments are calculated from the histogram. For constant bin width $\Delta x_{k}=\Delta x$ these two expressions can be simplified using $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}{\Big (}x_{k}-m_{1}^{(h)}{\Big )}^{n}h(x_{k})$

The second approach from Choi and Sweetman is an analytical methodology to combine statistical moments from individual segments of a time-history such that the resulting overall moments are those of the complete time-history. This methodology could be used for parallel computation of statistical moments with subsequent combination of those moments, or for combination of statistical moments computed at sequential times.

If $Q$ sets of statistical moments are known: $(\gamma _{0,q},\mu _{q},\sigma _{q}^{2},\alpha _{3,q},\alpha _{4,q})\quad$ for $q=1,2,\ldots ,Q$, then each $\gamma _{n}$ can be expressed in terms of the equivalent $n$ raw moments:

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

where $\gamma _{0,q}$ is generally taken to be the duration of the $q^{th}$ time-history, or the number of points if $\Delta t$ is constant.

The benefit of expressing the statistical moments in terms of $\gamma$ is that the $Q$ sets can be combined by addition, and there is no upper limit on the value of $Q$.

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

where the subscript $_{c}$ represents the concatenated time-history or combined $\gamma$. These combined values of $\gamma$ can then be inversely transformed into raw moments representing the complete concatenated time-history

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

Known relationships between the raw moments ($m_{n}$) and the central moments ($\theta _{n}=\operatorname {E} [(x-\mu )^{n}])$) are then used to compute the central moments of the concatenated time-history. Finally, the statistical moments of the concatenated history are computed from the central moments:

$\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 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}}.$

The two-pass algorithm may be written as:

def two_pass_covariance(data1, data2):
    n = len(data1)
    mean1 = sum(data1) / n
    mean2 = sum(data2) / n

    covariance = 0
    for i1, i2 in zip(data1, data2):
        a = i1 - mean1
        b = i2 - mean2
        covariance += a * b / n
    return covariance

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

Online

A stable one-pass algorithm exists, similar to the online algorithm for computing the variance, that computes 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}}$

The apparent asymmetry in that last equation is due to the fact that ${\textstyle (x_{n}-{\bar {x}}_{n})={\frac {n-1}{n}}(x_{n}-{\bar {x}}_{n-1})}$, so both update terms are equal to ${\textstyle {\frac {n-1}{n}}(x_{n}-{\bar {x}}_{n-1})(y_{n}-{\bar {y}}_{n-1})}$. Even greater accuracy can be achieved by first computing the means, then using the stable one-pass algorithm on the residuals.

Thus the covariance can be computed as

${\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}}$

def online_covariance(data1, data2):
    meanx = meany = C = n = 0
    for x, y in zip(data1, data2):
        n += 1
        dx = x - meanx
        meanx += dx / n
        meany += (y - meany) / n
        C += dx * (y - meany)

    population_covar = C / n
    # Bessel's correction for sample variance
    sample_covar = C / (n - 1)

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

A small modification can also be made to compute the weighted covariance:

def online_weighted_covariance(data1, data2, data3):
    meanx = meany = 0
    wsum = wsum2 = 0
    C = 0
    for x, y, w in zip(data1, data2, data3):
        wsum += w
        wsum2 += w * w
        dx = x - meanx
        meanx += (w / wsum) * dx
        meany += (w / wsum) * (y - meany)
        C += w * dx * (y - meany)

    population_covar = C / wsum
    # Bessel's correction for sample variance
    # Frequency weights
    sample_frequency_covar = C / (wsum - 1)
    # Reliability weights
    sample_reliability_covar = C / (wsum - wsum2 / wsum)

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

Likewise, there is a formula for combining the covariances of two sets that can be used to parallelize the computation:

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

Weighted batched version

A version of the weighted online algorithm that does batched updated also exists: let $w_{1},\dots w_{N}$ denote the weights, and write

${\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}}$

The covariance can then be computed as

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

Voir aussi

• Algorithme de sommation de Kahan
• Écart-type
• Méthode de Yamartino

References

External links

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