Algorithmes de calcul de la variance - Algorithms for calculating variance

Les algorithmes de calcul de la variance jouent un rôle majeur dans les statistiques computationnelles . Une difficulté clé dans la conception de bons algorithmes pour ce problème est que les formules de la variance peuvent impliquer des sommes de carrés, ce qui peut conduire à une instabilité numérique ainsi qu'à un débordement arithmétique lorsqu'il s'agit de grandes valeurs.

Algorithme naïf

Une formule pour calculer la variance d'une population entière de taille N est :

${\displaystyle \sigma ^{2}={\overline {(x^{2})}}-{\bar {x}}^{2}=\displaystyle {\frac {\sum _{i=1} ^{N}x_{i}^{2}-(\sum _{i=1}^{N}x_{i})^{2}/N}{N}}.\!}$

En utilisant la correction de Bessel pour calculer une estimation non biaisée de la variance de la population à partir d'un échantillon fini de n observations, la formule est :

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

Par conséquent, un algorithme naïf pour calculer la variance estimée est donné par ce qui suit :

Cet algorithme peut facilement être adapté pour calculer la variance d'une population finie : il suffit de diviser par N au lieu de n  − 1 sur la dernière ligne.

Étant donné que SumSq et (Sum×Sum)/ n peuvent être des nombres très similaires, l' annulation peut conduire à une précision du résultat bien inférieure à la précision inhérente de l' arithmétique à virgule flottante utilisée pour effectuer le calcul. Ainsi, cet algorithme ne doit pas être utilisé en pratique, et plusieurs algorithmes alternatifs, numériquement stables, ont été proposés. Ceci est particulièrement mauvais si l'écart type est faible par rapport à la moyenne.

Calcul de données décalées

La variance est invariante par rapport aux changements d'un paramètre d'emplacement , une propriété qui peut être utilisée pour éviter l'annulation catastrophique dans cette formule.

${\displaystyle \operatorname {Var} (XK)=\operatorname {Var} (X).}$

avec n'importe quelle constante, ce qui conduit à la nouvelle formule ${\style d'affichage K}$

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

plus le résultat est proche de la valeur moyenne, plus le résultat sera précis, mais le simple fait de choisir une valeur dans la plage des échantillons garantira la stabilité souhaitée. Si les valeurs sont petites, il n'y a pas de problèmes avec la somme de ses carrés, au contraire, si elles sont grandes, cela signifie nécessairement que la variance est également grande. Dans tous les cas, le deuxième terme de la formule est toujours plus petit que le premier, donc aucune annulation ne peut se produire. ${\style d'affichage K}$${\style d'affichage (x_{i}-K)}$

Si seul le premier échantillon est pris car l'algorithme peut être écrit en langage de programmation Python comme ${\style d'affichage K}$

def shifted_data_variance(data):
    if len(data) < 2:
        return 0.0
    K = data[0]
    n = Ex = Ex2 = 0.0
    for x in data:
        n = n + 1
        Ex += x - K
        Ex2 += (x - K) * (x - K)
    variance = (Ex2 - (Ex * Ex) / n) / (n - 1)
    # use n instead of (n-1) if want to compute the exact variance of the given data
    # use (n-1) if data are samples of a larger population
    return variance

Cette formule facilite également le calcul incrémental qui peut être exprimé sous la forme

K = n = Ex = Ex2 = 0.0

def add_variable(x):
    global K, n, Ex, Ex2
    if n == 0:
        K = x
    n += 1
    Ex += x - K
    Ex2 += (x - K) * (x - K)

def remove_variable(x):
    global K, n, Ex, Ex2
    n -= 1
    Ex -= x - K
    Ex2 -= (x - K) * (x - K)

def get_mean():
    global K, n, Ex
    return K + Ex / n

def get_variance():
    global n, Ex, Ex2
    return (Ex2 - (Ex * Ex) / n) / (n - 1)

Algorithme en deux passes

Une approche alternative, utilisant une formule différente pour la variance, calcule d'abord la moyenne de l'échantillon,

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

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

${\displaystyle {\text{sample variance}}=s^{2}=\displaystyle {\frac {\sum _{i=1}^{n}(x_{i}-{\bar {x}}) ^{2}}{n-1}},\!}$

où s est l'écart type. Ceci est donné par le code suivant :

def two_pass_variance(data):
    n = sum1 = sum2 = 0

    for x in data:
        n += 1
        sum1 += x

    mean = sum1 / n

    for x in data:
        sum2 += (x - mean) * (x - mean)

    variance = sum2 / (n - 1)
    return variance

Cet algorithme est numériquement stable si n est petit. Cependant, les résultats de ces deux algorithmes simples ("naïfs" et "à deux passes") peuvent dépendre excessivement de l'ordre des données et peuvent donner de mauvais résultats pour des ensembles de données très volumineux en raison d'erreurs d'arrondi répétées dans l'accumulation des données. sommes. Des techniques telles que la sommation compensée peuvent être utilisées pour combattre cette erreur dans une certaine mesure.

L'algorithme en ligne de Welford

Il est souvent utile de pouvoir calculer la variance en un seul passage , en inspectant chaque valeur une seule fois ; par exemple, lorsque les données sont collectées sans suffisamment de stockage pour conserver toutes les valeurs, ou lorsque les coûts d'accès à la mémoire dominent ceux du calcul. Pour un tel algorithme en ligne , une relation de récurrence est requise entre les quantités à partir desquelles les statistiques requises peuvent être calculées de manière numériquement stable. ${\style d'affichage x_{i}}$

Les formules suivantes peuvent être utilisées pour mettre à jour la moyenne et la variance (estimée) de la séquence, pour un élément supplémentaire x _n . Ici, désigne la moyenne d' échantillon des n premiers échantillons , leur variance d'échantillon biaisée et leur variance d'échantillon non biaisée . ${\displaystyle {\overline {x}}_{n}={\frac {1}{n}}\sum _{i=1}^{n}x_{i}}$${\style d'affichage (x_{1},\dots ,x_{n})}$${\displaystyle \sigma _{n}^{2}={\frac {1}{n}}\sum _{i=1}^{n}\left(x_{i}-{\overline {x} }_{n}\droit)^{2}}$${\displaystyle s_{n}^{2}={\frac {1}{n-1}}\sum _{i=1}^{n}\left(x_{i}-{\overline {x} }_{n}\droit)^{2}}$

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

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

${\displaystyle 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 souffrent d'instabilité numérique, car elles soustraient à plusieurs reprises un petit nombre d'un grand nombre qui s'échelonne avec n . Une meilleure quantité pour la mise à jour est la somme des carrés des différences par rapport à la moyenne actuelle, , notée ici : ${\displaystyle \textstyle \sum _{i=1}^{n}(x_{i}-{\bar {x}}_{n})^{2}}$${\style d'affichage M_{2,n}}$

${\displaystyle {\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 il a été minutieusement analysé. Il est également courant de désigner et . ${\displaystyle M_{k}={\bar {x}}_{k}}$${\style d'affichage S_{k}=M_{2,k}}$

Un exemple d'implémentation Python pour l'algorithme de Welford est donné ci-dessous.

# For a new value newValue, compute the new count, new mean, the new M2.
# mean accumulates the mean of the entire dataset
# M2 aggregates the squared distance from the mean
# count aggregates the number of samples seen so far
def update(existingAggregate, newValue):
    (count, mean, M2) = existingAggregate
    count += 1
    delta = newValue - mean
    mean += delta / count
    delta2 = newValue - mean
    M2 += delta * delta2
    return (count, mean, M2)

# Retrieve the mean, variance and sample variance from an aggregate
def finalize(existingAggregate):
    (count, mean, M2) = existingAggregate
    if count < 2:
        return float("nan")
    else:
        (mean, variance, sampleVariance) = (mean, M2 / count, M2 / (count - 1))
        return (mean, variance, sampleVariance)

Cet algorithme est beaucoup moins sujet à une perte de précision due à une annulation catastrophique , mais pourrait ne pas être aussi efficace en raison de l'opération de division à l'intérieur de la boucle. Pour un algorithme à deux passes particulièrement robuste pour le calcul de la variance, on peut d'abord calculer et soustraire une estimation de la moyenne, puis utiliser cet algorithme sur les résidus.

L' algorithme parallèle ci-dessous illustre comment fusionner plusieurs ensembles de statistiques calculées en ligne.

Algorithme incrémental pondéré

L'algorithme peut être étendu pour gérer des poids d'échantillon inégaux, en remplaçant le simple compteur n par la somme des poids vus jusqu'à présent. West (1979) propose cet algorithme incrémental :

def weighted_incremental_variance(data_weight_pairs):
    w_sum = w_sum2 = mean = S = 0

    for x, w in data_weight_pairs:
        w_sum = w_sum + w
        w_sum2 = w_sum2 + w * w
        mean_old = mean
        mean = mean_old + (w / w_sum) * (x - mean_old)
        S = S + w * (x - mean_old) * (x - mean)

    population_variance = S / w_sum
    # Bessel's correction for weighted samples
    # Frequency weights
    sample_frequency_variance = S / (w_sum - 1)
    # Reliability weights
    sample_reliability_variance = S / (w_sum - w_sum2 / w_sum)

Algorithme parallèle

Chan et al. notez que l'algorithme en ligne de Welford détaillé ci-dessus est un cas particulier d'algorithme qui fonctionne pour combiner des ensembles arbitraires et : ${\style d'affichage A}$${\style d'affichage B}$

${\displaystyle {\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} }}\\\fin{aligné}}}$.

Cela peut être utile lorsque, par exemple, plusieurs unités de traitement peuvent être affectées à des parties discrètes de l'entrée.

La méthode de Chan pour estimer la moyenne est numériquement instable lorsque et les deux sont grands, car l'erreur numérique n'est pas réduite comme elle l'est dans le cas. Dans de tels cas, préférez . ${\style d'affichage n_{A}\environ n_{B}}$${\displaystyle \delta ={\bar {x}}_{B}-{\bar {x}}_{A}}$${\style d'affichage n_{B}=1}$${\displaystyle {\bar {x}}_{AB}={\frac {n_{A}{\bar {x}}_{A}+n_{B}{\bar {x}}_{B} }{n_{AB}}}}$

def parallel_variance(n_a, avg_a, M2_a, n_b, avg_b, M2_b):
    n = n_a + n_b
    delta = avg_b - avg_a
    M2 = M2_a + M2_b + delta ** 2 * n_a * n_b / n
    var_ab = M2 / (n - 1)
    return var_ab

Cela peut être généralisé pour permettre la parallélisation avec AVX , avec les GPU , et les clusters d'ordinateurs , et à la covariance.

Exemple

Supposons que toutes les opérations à virgule flottante utilisent l' arithmétique double précision standard IEEE 754 . Considérons l'échantillon (4, 7, 13, 16) d'une population infinie. Sur la base de cet échantillon, la moyenne estimée de la population est de 10 et l'estimation sans biais de la variance de la population est de 30. L'algorithme naïf et l'algorithme à deux passes calculent correctement ces valeurs.

Considérons ensuite l'échantillon ( 10 ⁸  + 4 , 10 ⁸  + 7 , 10 ⁸  + 13 , 10 ⁸  + 16 ), qui donne lieu à la même variance estimée que le premier échantillon. L'algorithme à deux passes calcule correctement cette estimation de la variance, mais l'algorithme naïf renvoie 29.3333333333333332 au lieu de 30.

Bien que cette perte de précision puisse être tolérable et considérée comme un défaut mineur de l'algorithme naïf, l'augmentation supplémentaire du décalage rend l'erreur catastrophique. Considérons l'échantillon ( 10 ⁹  + 4 , 10 ⁹  + 7 , 10 ⁹  + 13 , 10 ⁹  + 16 ). Encore une fois, la variance de population estimée de 30 est calculée correctement par l'algorithme à deux passes, mais l'algorithme naïf la calcule maintenant comme −170,666666666666666. Ceci est un problème sérieux avec l'algorithme naïf et est dû à l' annulation catastrophique de la soustraction de deux nombres similaires à l'étape finale de l'algorithme.

Statistiques d'ordre supérieur

Terriberry étend les formules de Chan au calcul des troisième et quatrième moments centraux , nécessaires par exemple pour estimer l' asymétrie et l' aplatissement :

${\displaystyle {\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}\gauche(n_{A}^{2}-n_{A}n_{B}+n_{B}^{2}\droit)}{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{aligné}}}$

Voici à nouveau les sommes des puissances des différences par rapport à la moyenne , donnant ${\style d'affichage M_{k}}$${\displaystyle \sum (x-{\overline {x}})^{k}}$

${\displaystyle {\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émentiel (c'est-à-dire ), cela se simplifie en : ${\style d'affichage B=\{x\}}$

${\displaystyle {\begin{aligned}\delta &=xm\\[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 , une seule opération de division est nécessaire et les statistiques d'ordre supérieur peuvent ainsi être calculées pour un coût incrémental faible. ${\displaystyle \delta /n}$

Un exemple de l'algorithme en ligne pour l'aplatissement implémenté comme décrit est :

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 * delta_n
        term1 = delta * delta_n * n1
        mean = mean + delta_n
        M4 = M4 + term1 * delta_n2 * (n*n - 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 * M2) - 3
    return kurtosis

Pébaÿ étend encore ces résultats aux moments centraux d' ordre arbitraire , pour les cas incrémental et par paires, et par la suite Pébaÿ et al. pour les moments pondérés et composés. On peut aussi y trouver des formules similaires pour la covariance .

Choi et Sweetman proposent deux méthodes alternatives pour calculer l'asymétrie et l'aplatissement, chacune pouvant économiser des exigences substantielles en matière de mémoire informatique et de temps CPU dans certaines applications. La première approche consiste à calculer les moments statistiques en séparant les données en bacs, puis en calculant les moments à partir de la géométrie de l'histogramme résultant, qui devient effectivement un algorithme à une passe pour les moments plus élevés. Un avantage est que les calculs de moment statistique peuvent être effectués avec une précision arbitraire de sorte que les calculs peuvent être ajustés à la précision, par exemple, du format de stockage de données ou du matériel de mesure d'origine. Un histogramme relatif d'une variable aléatoire peut être construit de la manière conventionnelle : la plage de valeurs potentielles est divisée en cases et le nombre d'occurrences dans chaque case est compté et tracé de telle sorte que l'aire de chaque rectangle soit égale à la partie des valeurs de l'échantillon dans ce bac :

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

où et représentent la fréquence et la fréquence relative à bin et est la surface totale de l'histogramme. Après cette normalisation, les moments bruts et moments centraux de peuvent être calculés à partir de l'histogramme relatif : ${\style d'affichage h(x_{k})}$${\style d'affichage H(x_{k})}$${\style d'affichage x_{k}}$${\displaystyle A=\sum _{k=1}^{K}h(x_{k})\,\Delta x_{k}}$${\style d'affichage n}$${\style d'affichage x(t)}$

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

${\displaystyle \theta _{n}^{(h)}=\sum _{k=1}^{K}{\Big (}x_{k}-m_{1}^{(h)}{\ Grand )}^{n}\,H(x_{k})\,\Delta x_{k}={\frac {1}{A}}\sum _{k=1}^{K}{\Grand (}x_{k}-m_{1}^{(h)}{\Grand )}^{n}h(x_{k})\,\Delta x_{k}}$

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

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

${\displaystyle \theta _{n}^{(h)}={\frac {1}{I}}\sum _{k=1}^{K}{\Big (}x_{k}-m_{ 1}^{(h)}{\Grand )}^{n}h(x_{k})}$

La deuxième approche de Choi et Sweetman est une méthodologie analytique pour combiner des moments statistiques de segments individuels d'une histoire temporelle de telle sorte que les moments globaux résultants soient ceux de l'histoire temporelle complète. Cette méthodologie pourrait être utilisée pour le calcul parallèle de moments statistiques avec une combinaison ultérieure de ces moments, ou pour une combinaison de moments statistiques calculés à des moments séquentiels.

Si des ensembles de moments statistiques sont connus : pour , alors chacun peut être exprimé en termes de moments bruts équivalents : ${\style d'affichage Q}$${\displaystyle (\gamma _{0,q},\mu _{q},\sigma _{q}^{2},\alpha _{3,q},\alpha _{4,q})\ quad }$${\displaystyle q=1,2,\ldots ,Q}$${\displaystyle \gamma _{n}}$${\style d'affichage n}$

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

où est généralement considéré comme étant la durée de l' historique temporel, ou le nombre de points si est constant. ${\displaystyle \gamma _{0,q}}$${\displaystyle q^{th}}$${\style d'affichage \Delta t}$

L'avantage d'exprimer les moments statistiques en termes de est que les ensembles peuvent être combinés par addition et qu'il n'y a pas de limite supérieure à la valeur de . ${\style d'affichage \gamma }$${\style d'affichage Q}$${\style d'affichage Q}$

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

où l'indice représente l'historique temporel concaténé ou combiné . Ces valeurs combinées de peuvent ensuite être inversement transformées en moments bruts représentant l'histoire temporelle concaténée complète ${\style d'affichage _{c}}$${\style d'affichage \gamma }$${\style d'affichage \gamma }$

${\displaystyle 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 ( ) et les moments centraux ( ) sont ensuite utilisées pour calculer les moments centraux de l'histoire temporelle concaténée. Enfin, les moments statistiques de l'histoire concaténée sont calculés à partir des moments centraux : ${\style d'affichage m_{n}}$${\displaystyle \theta _{n}=\operatorname {E} [(x-\mu )^{n}])}$

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

Covariance

Des algorithmes très similaires peuvent être utilisés pour calculer la covariance .

Algorithme naïf

L'algorithme naïf est :

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

Pour l'algorithme ci-dessus, on pourrait utiliser le code Python suivant :

def naive_covariance(data1, data2):
    n = len(data1)
    sum12 = 0
    sum1 = sum(data1)
    sum2 = sum(data2)

    for i1, i2 in zip(data1, data2):
        sum12 += i1 * i2

    covariance = (sum12 - sum1 * sum2 / n) / n
    return covariance

Avec estimation de la moyenne

En ce qui concerne la variance, la covariance de deux variables aléatoires est également invariante par décalage, donc étant donné deux valeurs constantes et cela peut s'écrire : ${\style d'affichage k_{x}}$${\displaystyle k_{y},}$

${\displaystyle \operatorname {Cov} (X,Y)=\operatorname {Cov} (X-k_{x},Y-k_{y})=\displaystyle {\frac {\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 encore une fois, le choix d'une valeur à l'intérieur de la plage de valeurs stabilisera la formule contre une annulation catastrophique et la rendra plus robuste contre les grosses sommes. En prenant la première valeur de chaque ensemble de données, l'algorithme peut être écrit comme :

def shifted_data_covariance(data_x, data_y):
    n = len(data_x)
    if n < 2:
        return 0
    kx = data_x[0]
    ky = data_y[0]
    Ex = Ey = Exy = 0
    for ix, iy in zip(data_x, data_y):
        Ex += ix - kx
        Ey += iy - ky
        Exy += (ix - kx) * (iy - ky)
    return (Exy - Ex * Ey / n) / n

Deux passes

L'algorithme à deux passes calcule d'abord les moyennes de l'échantillon, puis la covariance :

${\displaystyle {\bar {x}}=\displaystyle \sum _{i=1}^{n}x_{i}/n}$

${\displaystyle {\bar {y}}=\displaystyle \sum _{i=1}^{n}y_{i}/n}$

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

L'algorithme à deux passes peut s'écrire :

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

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 devraient être nulles, mais la deuxième passe compense toute petite erreur. ${\displaystyle \textstyle \sum x_{i}}$${\displaystyle \textstyle \sum y_{i}}$

En ligne

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

${\displaystyle {\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 cette dernière équation est due au fait que , donc les deux termes de mise à jour sont égaux à . Une précision encore plus grande peut être obtenue en calculant d'abord les moyennes, puis en utilisant l'algorithme stable à un passage sur les résidus. ${\displaystyle \textstyle (x_{n}-{\bar {x}}_{n})={\frac {n-1}{n}}(x_{n}-{\bar {x}}_ {n-1})}$${\displaystyle \textstyle {\frac {n-1}{n}}(x_{n}-{\bar {x}}_{n-1})(y_{n}-{\bar {y}} _{n-1})}$

Ainsi, la covariance peut être calculée comme

${\displaystyle {\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)

Une petite modification peut également être apportée pour calculer la covariance pondérée :

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)

De même, il existe une formule pour combiner les covariances de deux ensembles qui peut être utilisée pour paralléliser le calcul :

${\displaystyle 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 par lots pondérée

Une version de l'algorithme en ligne pondéré qui fait la mise à jour par lots existe également : let dénoter les poids, et écrire ${\displaystyle w_{1},\dots w_{N}}$

${\displaystyle {\begin{alignedat}{2}{\bar {x}}_{n+k}&={\bar {x}}_{n}&\,+\,&{\frac {\ somme _{i=n+1}^{n+k}w_{i}(x_{i}-{\bar {x}}_{n})}{\somme _{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 peut alors être calculée comme

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

Voir également

• Algorithme de sommation de Kahan
• Écarts au carré par rapport à la moyenne
• Méthode Yamartino

Les références

Liens externes

• Weisstein, Eric W. "Calcul de la variance d'échantillon" . MathWorld .