Les algorithmes de calcul de variance - Algorithms for calculating variance


Un article de Wikipédia, l'encyclopédie libre

Les algorithmes de calcul de variance jouent un rôle majeur dans les statistiques de calcul . Une difficulté clé dans la conception des bons algorithmes pour ce problème est que les formules de la variance peuvent comporter des sommes des carrés, ce qui peut conduire à l' instabilité numérique , ainsi que de dépassement arithmétique lorsqu'ils traitent avec de grandes valeurs.

algorithme naïf

Une formule de calcul de la variance d'un ensemble de population de taille N est:

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

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

  • Let n ← 0, Somme ← 0, SUMSQ ← 0
  • Pour chaque donnée x :
    • nn + 1
    • Somme ← Somme + x
    • SUMSQ ← SUMSQ + x x x
  • Var = (SUMSQ - (Sum Somme x) / n) / (n - 1)

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.

Parce que SUMSQ et (Somme × Somme) / n peuvent être des nombres très similaires, l' annulation peut conduire à la précision du résultat soit beaucoup moins que la précision inhérente de l' arithmétique en virgule flottante utilisé pour effectuer le calcul. Ainsi , cet algorithme ne doit pas être utilisé dans la pratique, et plusieurs autre, numériquement stables, les algorithmes ont été proposés. Ceci est particulièrement mauvais si l'écart type est faible par rapport à la moyenne. Cependant, l'algorithme peut être améliorée en adoptant la méthode de la moyenne supposée .

Computing données décalées

Nous pouvons utiliser une propriété de la variance pour éviter l'annulation catastrophique dans cette formule, à savoir la variance est invariante par rapport aux changements d'un paramètre de localisation

avec une constante quelconque, ce qui conduit à la nouvelle formule

le plus proche est à la valeur moyenne plus le résultat sera précis, mais juste le choix d' une valeur à l' intérieur des échantillons vont garantira la stabilité souhaitée. Si les valeurs sont faibles alors il n'y a aucun problème avec la somme de ses carrés, au contraire, si elles sont grandes , cela signifie nécessairement que la variance est importante aussi bien. Dans tous les cas , le second terme de la formule est toujours plus petit que le premier donc aucune annulation peut se produire.

Si nous prenons juste le premier échantillon que l'algorithme peut être écrit en langage de programmation Python comme

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 ainsi le calcul différentiel, qui peut être exprimé en

K = n = Ex = Ex2 = 0.0

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

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

def get_meanvalue():
    return K + Ex / n

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

Deux passes algorithme

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

,

et calcule ensuite la somme des carrés des différences de la moyenne,

,

où s est l'écart-type. Ceci est donné par le pseudo-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 faible. Cependant, les résultats de ces deux algorithmes simples ( « Naïf » et « deux passes ») peuvent compter démesurément sur l'ordre des données et peuvent donner de mauvais résultats pour les ensembles de données très volumineux en raison d' une erreur roundoff répétée dans l'accumulation du sommes. Des techniques telles que la somme compensée peuvent être utilisés pour lutter contre cette erreur dans une certaine mesure.

algorithme en ligne de Welford

Il est souvent utile de pouvoir calculer la variance en une seule passe, inspectant chaque valeur une seule fois; par exemple, lorsque les données sont collectées sans stockage suffisant pour garder toutes les valeurs, ou lorsque les coûts d'accès à la mémoire dominent ceux de calcul. Pour un tel algorithme en ligne , une relation de récurrence est nécessaire entre les quantités dont les statistiques nécessaires peuvent être calculées de façon stable numériquement.

Les formules suivantes peuvent être utilisées pour mettre à jour la moyenne de variance et (estimation) de la séquence, pour un élément supplémentaire x n . Ici, x n désigne la moyenne de l' échantillon des premiers n échantillons ( x 1 , ..., x n ), s 2 n leur variance de l' échantillon et σ 2 n leur variance de la population.

Ces formules souffrent d' une instabilité numérique, comme ils soustraient à plusieurs reprises un petit nombre d'un grand nombre qui évolue avec n . Une meilleure quantité pour la mise à jour est la somme des carrés des différences par rapport à la moyenne actuelle, , ici notée :

Cet algorithme a été trouvé par Welford, et il a été analysé en profondeur. Il est également courant pour désigner et .

Un exemple de mise en œuvre Python pour l'algorithme de Welford est donnée 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
    (mean, variance, sampleVariance) = (mean, M2/count, M2/(count - 1)) 
    if count < 2:
        return float('nan')
    else:
        return (mean, variance, sampleVariance)

Cet algorithme est beaucoup moins sujette à la perte de précision en raison de l' 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 passe particulièrement robuste pour calculer la variance, on peut d' abord calculer et soustraire l'estimation de la moyenne, et ensuite utiliser cet algorithme sur les résidus.

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

Pondérée algorithme incrémental

L'algorithme peut être étendu pour gérer le poids des échantillons inégales, en remplaçant le simple compteur n avec la somme des poids observés jusqu'à présent. Ouest (1979) suggère cet algorithme incrémental :

def weighted_incremental_variance(dataWeightPairs):
    wSum = wSum2 = mean = S = 0

    for x, w in dataWeightPairs:  # Alternatively "for x, w in zip(data, weights):"
        wSum = wSum + w
        wSum2 = wSum2 + w*w
        meanOld = mean
        mean = meanOld + (w / wSum) * (x - meanOld)
        S = S + w * (x - meanOld) * (x - mean)

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

algorithme parallèle

Chan et al. noter que l'algorithme ci - dessus « en ligne » est un cas particulier d'un algorithme qui fonctionne pour une partition de l'échantillon en ensembles , :

.

Ceci peut être utile lorsque, par exemple, des unités de traitement multiples peuvent être affectés à 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 grandes, parce que l'erreur numérique dans n'est pas réduite de la manière qu'il est dans le cas. Dans ce cas, préférez .

def parallel_variance(avg_a, count_a, var_a, avg_b, count_b, var_b):
    delta = avg_b - avg_a
    m_a = var_a * (count_a - 1)
    m_b = var_b * (count_b - 1)
    M2 = m_a + m_b + delta ** 2 * count_a * count_b / (count_a + count_b)
    return M2 / (count_a + count_b - 1)

Exemple

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

Considérons maintenant l'échantillon ( 10 8  + 4 , 10 8  + 7 , 10 8  + 13 , 10 8  + 16 ), ce qui donne lieu à la même variance estimée comme étant le premier échantillon. L'algorithme calcule deux passes cette estimation de la variance correctement, mais l'algorithme naïf retourne 29,333333333333332 au lieu de 30.

Bien que cette perte de précision peut être tolérable et considéré comme un défaut mineur de l'algorithme naïf, ce qui augmente encore le décalage fait l'erreur catastrophique. Considérons l'échantillon ( 10 9  + 4 , 10 9  + 7 , 10 9  + 13 , 10 9  + 16 ). Encore une fois la variance de la population estimée de 30 est calculée correctement par l'algorithme à deux passes, mais l'algorithme naïf , il calcule maintenant -170,66666666666666. Ceci est un problème sérieux avec l' algorithme naïf et est due à l' annulation catastrophique dans la soustraction de deux nombres similaires à l'étape finale de l'algorithme.

Statistiques d'ordre supérieur

Terriberry étend les formules de calcul Chan pour les troisième et quatrième moments centraux , nécessaire par exemple lors de l' estimation d' asymétrie et de kurtosis :

Ici , le sont à nouveau les sommes des puissances des différences par rapport à la moyenne , ce qui donne

dissymétrie:
kurtosis:

Pour le cas incrémental (c. -à ), ce qui simplifie à:

En préservant la valeur , une seule opération de division est nécessaire et les statistiques d'ordre supérieur peut ainsi être calculé pour peu de frais supplémentaires.

Un exemple de l'algorithme en ligne pour kurtosis mis en œuvre comme décrit est le suivant:

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

    kurtosis = (n*M4) / (M2*M2) - 3
    return kurtosis

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

Choi et Sweetman offrent deux méthodes alternatives pour calculer la skewness et kurtosis, chacun peut enregistrer les besoins en mémoire d'ordinateur et substantielle temps CPU dans certaines applications. La première approche consiste à calculer les moments statistiques en séparant les données dans les bacs, puis le calcul des moments de la géométrie de l'histogramme obtenu, qui devient effectivement un algorithme en une seule passe pour des moments plus élevés. L' un des avantages est que les calculs de moment statistique peuvent être réalisées avec une précision arbitraire de telle sorte que les calculs peuvent être accordés à la précision, par exemple, le format de stockage de données ou le matériel de mesure d' origine. Un histogramme relatif d'une variable aléatoire peut être construit de façon classique: la plage de valeurs possibles est divisée en cases et le nombre d'occurrences de chaque bac sont comptées et tracé de telle sorte que la surface de chaque rectangle est égale à la partie des valeurs d'échantillon dans ce bac:

où et représentent la fréquence et la fréquence relative au casier et est la surface totale de l'histogramme. Après cette normalisation, les moments premières et des moments centraux de peuvent être calculées à partir de l'histogramme relatif:

où l'indice indique les moments sont calculés à partir de l'histogramme. Pour la largeur du bac constant peut être simplifié ces deux expressions en utilisant :

La deuxième approche de Choi et Sweetman est une méthodologie analytique de combiner des moments statistiques de segments individuels d'un temps historique telle que les moments globaux résultant sont ceux du temps complet de l'histoire. Cette méthode pourrait être utilisée pour le calcul parallèle des moments statistiques avec la combinaison ultérieure de ces moments, ou pour la combinaison des moments statistiques calculés à des moments successifs.

Si des ensembles de moments statistiques sont connus: pour , chacun peut être exprimé en termes des équivalents moments premières:

où est généralement considéré comme la durée du temps historique, ou le nombre de points si est constant.

L'avantage d'exprimer les moments statistiques en termes de est que les ensembles peuvent être combinés par addition, et il n'y a pas de limite supérieure à la valeur de .

où l'indice représente la concaténés temps historique ou combiné . Ces valeurs combinées de peut alors être inversement transformées en moments premières représentant le temps complet concaténés histoire

Relations connues entre les moments premières ( ) et les moments centraux ( ) sont ensuite utilisés pour calculer les moments centraux du concaténés temps historique. Enfin, les moments statistiques de l'histoire concaténés sont calculées à partir des moments centraux:

covariance

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

algorithme naïf

L'algorithme est naïf:

Pour l'algorithme ci-dessus, on peut 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 invariante déplace aussi, donc donné deux valeurs constantes et il peut être écrit:

et en choisissant à nouveau une valeur dans la plage des valeurs stabilisera la formule contre l'annulation catastrophique, ainsi que le rendre plus robuste contre de grosses sommes. En prenant la première valeur de chaque ensemble de données, l'algorithme peut être écrit:

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

Deux passes

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

L'algorithme en deux passes peut être écrit:

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é un peu plus précise effectue l'algorithme complet naïf sur les résidus. Les sommes finales et devrait être égale à zéro, mais la deuxième passe compense toute petite erreur.

en ligne

Un algorithme d' une passe stable existe, semblable à l'algorithme en ligne pour calculer la variance, qui calcule co-moment :

L'asymétrie apparente dans cette dernière équation est due au fait que , si les deux termes de mise à jour sont égales . Même une plus grande précision peut être obtenue en calculant d' abord les moyens, puis en utilisant la stable en un seul passage algorithme sur les résidus.

On peut donc calculer la covariance comme

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)

Nous pouvons également faire une petite modification 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 des deux ensembles qui peuvent être utilisés pour paralléliser le calcul:

Version pondérée en lot

Une version de l'algorithme pondéré en ligne qui n'existe aussi mis à jour par lots: laisser désigner les poids, nous pouvons écrire

La covariance peut alors être calculé comme

Voir également

Références

  1. ^ Un b Bo Einarsson (1 Août 2005). Précision et fiabilité en informatique scientifique . SIAM. p. 47. ISBN  978-0-89871-584-2 . Récupéré 17 Février 2013 .
  2. ^ Un b T.F.Chan, GH Golub et RJ LeVeque (1983). « » Les algorithmes de calcul de la variance de l' échantillon: analyse et recommandations », le Statisticien américain, 37" (PDF) : 242-247.
  3. ^ A b Schubert, Erich; Gertz, Michael (09/07/2018). « Calcul parallèle Numériquement stable de la variance (co) » . ACM: 10. doi : 10,1145 / 3.221.269,3223036 . ISBN  9781450365055 .
  4. ^ Higham, Nicolas (2002). La précision et la stabilité des algorithmes numériques (2) (éd problème 1.10) . SIAM.
  5. ^ BP Welford (1962). « Note sur une méthode de calcul des sommes corrigées des carrés et des produits » . Technométrie 4 (3): 419-420.
  6. ^ Donald E. Knuth (1998). L'art de la programmation informatique , volume 2: Seminumerical algorithmes ., 3e EDN, p. 232. Boston: Addison-Wesley.
  7. ^ Chan, Tony F .; Golub, Gene H .; LeVeque, J. Randall (1983). Les algorithmes de calculla variance deéchantillon: Analyse et recommandations. Le Statisticien américain 37, 242-247. https://www.jstor.org/stable/2683386
  8. ^ Ling, Robert F. (1974). Comparaison de plusieurs algorithmes pourmoyenscalcul de Echantillon et Variances. Journal de l'American Statistical Association, Vol. 69, n ° 348, 859-866. doi : 10,2307 / 2286154
  9. ^ http://www.johndcook.com/standard_deviation.html
  10. ^ DHD Ouest (1979). Communications of the ACM , 22, 9, 532-535: Mise à jour des estimations moyenne et la variance: un procédé amélioré
  11. ^ Chan, Tony F. ; Golub, Gene H. ; LeVeque, J. Randall (1979), "Mise à jour des formules et un algorithme pour le calcul des échantillons par paires Variances." (PDF) , Rapport technique STAN-CS-79-773 , Département d'informatique, Université de Stanford.
  12. ^ Terriberry, Timothy B. (2007), Informatique ordre supérieur Moments en ligne
  13. ^ PEBAY, Philippe (2008), "Formules pour Robuste, One-Pass Calcul parallèle de Covariances et arbitraires ordre Moments statistiques" (PDF) , Rapport technique SAND2008-6212 , Sandia National Laboratories
  14. ^ PEBAY, Philippe; Terriberry, Timothy; Kolla, Hemanth; Bennett, Janine (2016), « Numériquement stables, formules évolutives pour le calcul parallèle et en ligne d' ordre supérieur de multivariée moments centraux avec des poids arbitraires » , statistiques de calcul , Springer
  15. ^ A b Choi, Myoungkeun; Sweetman, Bert (2010), "Calcul efficace des moments statistiques pour Structural Health Monitoring" , Journal of Structural Health Monitoring , 9 (1)

Liens externes