Algoritmi pentru calcularea varianței

De la Wikipedia, enciclopedia liberă.
Salt la navigare Salt la căutare

Algoritmii pentru calcularea varianței joacă un rol foarte important în statisticile de calcul . O dificultate cheie în proiectarea unui algoritm bun pentru această problemă este că formulele pentru varianță pot include sume de pătrate, ceea ce poate duce la instabilitate numerică , precum și la deversare aritmetică atunci când se tratează valori mari.

Algoritm banal

O formulă pentru calcularea varianței unei populații întregi de dimensiunea N este:

Folosind corecția lui Bessel pentru a calcula o estimare fără părtinire a varianței populației dintr-un eșantion finit de n observații, formula este:

Prin urmare, un algoritm banal pentru calcularea estimării varianței este dat de:

 Def estimare Varianță ( x ):
   n = lungime ( x )
   sumă = 0
   sumQuad = 0
   pentru xi în x :
      sum = sum + xi
      sumQuad = sumQuad + xi ^ 2
   return ( sumQuad - sum ^ 2 / n ) / ( n - 1 )

Acest algoritm poate fi ușor adaptat pentru a calcula varianța unei populații finite: împărțiți doar cu N în loc de n - 1 în ultima linie.

Atâta timp cât Și pot fi numere foarte similare, ștergerea poate face ca precizia rezultatului să fie mult mai mică decât precizia intrinsecă a aritmeticii în virgulă mobilă utilizată pentru efectuarea calculului. Astfel, acest algoritm nu trebuie utilizat în practică. [1] [2] Acest lucru este deosebit de dezavantajos dacă deviația standard este mică în comparație cu media. Cu toate acestea, algoritmul poate fi îmbunătățit adoptând metoda medie ipotezată.

Calcul cu date traduse

Putem folosi o proprietate de varianță pentru a evita anularea numerică în această formulă, deoarece varianța este invariantă în ceea ce privește modificările modificărilor într-un parametru de locație

cu constant, ceea ce duce la noua formulă

și cu cât este mai aproape la valoarea medie, cu atât rezultatul va fi mai precis, dar doar alegerea unei valori în intervalul eșantionului va garanta stabilitatea dorită. Dacă valorile sunt mici, atunci nu există nicio problemă cu suma pătratelor lor; dimpotrivă, dacă sunt mari, aceasta înseamnă neapărat că și varianța este mare. În orice caz, al doilea termen din formulă este întotdeauna mai mic decât primul, deci nu poate avea loc nicio anulare. [2]

Dacă luăm doar prima mostră ca. , algoritmul poate fi scris în limbajul de programare Python după cum urmează

 def shifted_data_variance ( data ):
   dacă len ( date ) == 0 :
      returnează 0
   K = data [ 0 ]
   n = 0
   sum_ = 0
   sum_sqr = 0
   pentru data x în :
      n = n + 1
      sum_ + = x - K
      sum_sqr + = ( x - K ) * ( x - K )
   varianță = ( sum_sqr - ( sum_ * sum_ ) / n ) / ( n - 1 )
   # folosiți n în loc de (n-1) dacă doriți să calculați varianța exactă a datelor date
   # utilizați (n-1) dacă datele sunt eșantioane ale unei populații mai mari
   returnează varianța

această formulă a facilitat astfel calculul incremental, care poate fi exprimat după cum urmează

 K = 0
n = 0
ex = 0
ex2 = 0
 
def add_variable ( x ):
    dacă ( 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 ():
   returnează K + ex / n

def get_variance ():
    returnare ( ex2 - ( ex * ex ) / n ) / ( n - 1 )

Algoritm în doi pași

O abordare alternativă, folosind o formulă diferită pentru varianță, este să calculați mai întâi media eșantionului,

,

și apoi calculați suma pătratelor diferențelor față de medie,

,

unde s este abaterea standard. Acest lucru este dat de următorul cod sursă:

 def two_pass_variance ( data ):
    n = 0
    sum1 = 0
    sum2 = 0
    
    pentru data x în :
        n + = 1
        sum1 + = x
    
    medie = sum1 / n

    pentru data x în :
        sum2 + = ( x - medie ) * ( x - medie )
    
    varianță = sum2 / ( n - 1 )
    returnează varianța

Acest algoritm este stabil din punct de vedere numeric dacă n este mic. [1] [3] Cu toate acestea, rezultatele atât ale acestor algoritmi simpli („trivial”, cât și „cu două treceri”) pot depinde necomutativ de sortarea datelor și pot da rezultate slabe pentru seturi de date mari din cauza erorilor de rotunjiri repetate care se acumulează în sume. Tehnici precum suma compensată pot fi utilizate pentru a reduce această eroare.

Varianta compensată

Versiunea compensată sumar a algoritmului de mai sus este: [4]

 def Varianta compensată ( dată ):
    n = 0
    sum1 = 0
    pentru data x în :
        n + = 1
        sum1 + = x
    medie = sum1 / n
     
    sum2 = 0
    sum3 = 0
    pentru data x în :
        sum2 + = ( x - medie ) ** 2
        sum3 + = ( x - medie )
    varianță = ( sum2 - sum3 ** 2 / n ) / ( n - 1 )
    returnează varianța

Algoritm online

Este adesea util să puteți calcula varianța într-un singur pas, procesând fiecare valoare o singura data; de exemplu, atunci când datele sunt colectate fără spațiu suficient pentru a păstra toate valorile sau când costurile de accesare a memoriei sunt mai mari decât cele din calcul. Pentru un astfel de algoritm online, este necesară o relație de recurență între cantități din care statisticile necesare pot fi calculate într-un mod numeric stabil.

Următoarele formule pot fi utilizate pentru a actualiza media și varianța (estimată) a secvenței pentru un element suplimentar . Aici, x n reprezintă media eșantionului primelor n eșantioane ( x 1 , ..., x n ), s 2 n varianța eșantionului lor și σ 2 n varianța populației lor.

Aceste formule suferă de instabilitate numerică. O cantitate mai bună pentru actualizare este suma pătratelor diferențelor față de media curentă , aici indicat ca :

Un algoritm stabil din punct de vedere numeric pentru varianța eșantionului este prezentat mai jos. De asemenea, calculează media. Acest algoritm se datorează lui Knuth, [5] care îl citează pe Welford, [6] și a fost analizat cu atenție. [7] [8] Este, de asemenea, obișnuit să denotați Și . [9]

 def online_variance ( data ):
    n = 0
    medie = 0,0
    M2 = 0,0
     
    pentru data x în :
        n + = 1
        delta = x - medie
        medie + = delta / n
        M2 + = delta * ( x - medie )

    dacă n < 2 :
        return float ( 'nan' )
    altceva :
        retur M2 / ( n - 1 )

Acest algoritm este mult mai puțin predispus la pierderea preciziei datorită anulării numerice , dar este posibil să nu fie eficient datorită operației de divizare din buclă. Pentru un algoritm în doi pași deosebit de robust pentru calcularea varianței, se poate calcula și scădea mai întâi o estimare a mediei și apoi utiliza acest algoritm pe reziduuri.

Mai jos, algoritmul paralel ilustrează cum se îmbină mai multe seturi de statistici calculate online.

Algoritm incremental ponderat

Algoritmul poate fi extins pentru a gestiona greutăți eșantion inegale prin înlocuirea contorului simplu n cu suma greutăților văzute până acum. West (1979) [10] sugerează acest algoritm incremental:

 def ponderată_varianță_incrementală ( dataWeightPairs ):
    sumă = 0
    medie = 0
    M2 = 0

    pentru x , greutate în dataWeightPairs : # Alternativ "pentru x, greutate în zip (date, greutăți):"
        temp = greutate + greutate totală
        delta = x - medie
        R = delta * greutate / temp
        medie + = R
        M2 + = greutate * delta * R # Alternativ, „M2 = M2 + greutate * delta * (x - medie)”
        sumweight = temp

    varianță_n = M2 / greutate sumă
    varianță = varianță_n * len ( dataWeightPairs ) / ( len ( dataWeightPairs ) - 1 )

Algoritm paralel

Chan și alții. [4] Rețineți că algoritmul în linie prezentat mai sus este un caz special al unui algoritm care funcționează pentru fiecare partiție a eșantionului în seturi , :

.

Acest lucru poate fi util atunci când, de exemplu, mai multe unități de procesare pot fi atribuite unor bucăți de intrare discrete.

Metoda lui Chan de estimare a mediei este instabilă numeric atunci când și ambele sunt mari, deoarece eroarea numerică din nu este redus pentru a fi în cauză . În astfel de cazuri, este de preferat ca. .

Exemplu

Să presupunem că toate operațiunile în virgulă mobilă utilizează aritmetica standard IEEE 754 cu precizie dublă . Să luăm în considerare eșantionul (4, 7, 13, 16) dintr-o populație infinită. Pe baza acestui eșantion, media populației estimate este de 10, iar estimarea fără părtinire a varianței populației este de 30. Atât algoritmii „triviali”, cât și cei „în doi pași” calculează corect aceste valori. Acum, în schimb, să luăm în considerare eșantionul (10 8 + 4, 10 8 + 7, 10 8 + 13, 10 8 + 16), care dă naștere la aceeași varianță estimată ca pentru primul eșantion. Algoritmul „două treceri” calculează corect această varianță estimată, dar algoritmul „banal” returnează 29,333333333333332 în loc de 30. În timp ce această pierdere de precizie poate fi tolerabilă și considerată ca un defect minor al algoritmului „două treceri”, este ușor să găsiți date care relevă un defect grav în algoritmul „banal”: Să luăm eșantionul (10 9 + 4, 10 9 + 7, 10 9 + 13, 10 9 + 16). Din nou, varianța estimată a populației, de valoarea 30, este corect calculată de algoritmul „în doi pași”, dar algoritmul „banal” calculează acum −170,66666666666666. Aceasta este o problemă serioasă cu algoritmul „banal” și se datorează anulării numerice în scăderea a două numere similare în etapa finală a algoritmului.

Statistici de ordin superior

Terriberry [11] extinde formulele lui Chan la calcularea celui de-al treilea și al patrulea moment central, necesare de exemplu atunci când se estimează asimetria și curtoză :

Aici sunt din nou sumele puterilor diferențelor față de medie , având

asimetrie:
cutroza:

Pentru cazul incremental (de exemplu, ), aceasta se reduce la:

Păstrarea valorii , este necesară o singură operațiune de diviziune și statistici de ordin superior pot fi astfel calculate pentru un cost incremental mic.

Un exemplu de algoritm în linie pentru kurtoză implementat așa cum este descris este:

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

    pentru data x în :
        n1 = n
        n = n + 1
        delta = x - medie
        delta_n = delta / n
        delta_n2 = delta_n * delta_n
        termen1 = delta * delta_n * n1
        medie = medie + delta_n
        M4 = M4 + termen1 * delta_n2 * ( n * n - 3 * n + 3 ) + 6 * delta_n2 * M2 - 4 * delta_n * M3
        M3 = M3 + termen1 * delta_n * ( n - 2 ) - 3 * delta_n * M2
        M2 = M2 + termen1

    kurtosis = ( n * M4 ) / ( M2 * M2 ) - 3
    întoarce kurtosis

Pébay [12] extinde în continuare aceste rezultate pentru momente centrale de ordine arbitrară, pentru cazul incremental și pentru cazul cuplurilor. Formule similare pentru covarianță pot fi găsite.

Choi și Sweetman [13] oferă două metode alternative pentru calcularea asimetriei și a curtozei, fiecare dintre acestea putând economisi cerințele de memorie ale computerului și timpul de ocupare a procesorului în anumite aplicații. Prima abordare este de a calcula momentele statistice separând datele în pachete și apoi de a calcula momentele din geometria histogramei rezultată, care devine efectiv un algoritm cu un singur pas pentru momentele superioare. Un avantaj este că calculul momentelor statistice se poate face cu o precizie arbitrară, astfel încât calculele să poată fi aduse, de exemplu, la acuratețea formatului de stocare a datelor sau la capacitatea inițială a hardware-ului. O histogramă legată de o variabilă aleatorie poate fi construită în mod convențional: gama valorilor potențiale este împărțită în pachete și numărul de apariții din fiecare pachet este numărat și trasat astfel încât aria fiecărui dreptunghi să fie egală cu porțiunea de valori eșantion din acest pachet:

unde este Și reprezintă frecvența și frecvența relativă pentru pachet Și este aria totală a histogramei. După această normalizare, momentele de origine zero (momente brute) și momentele centrale ale poate fi calculat din histograma relevantă:

unde îl suprascrie indică faptul că momentele sunt calculate din histogramă. Pentru un pachet cu lățime constantă aceste două expresii pot fi simplificate folosind :

A doua abordare, de la Choi și Sweetman [13] , constă dintr-o metodologie analitică pentru combinarea momentelor statistice din segmente individuale ale istoriei timpului, astfel încât momentele globale rezultate să fie cele ale unei istorii complete a timpului. Această metodologie ar putea fi utilizată pentru calculul paralel al momentelor statistice cu combinația ulterioară a acelor momente sau pentru combinația momentelor statistice calculate în timpi secvențiali.

Dacă seturile de momente statistice sunt cunoscute: pentru , apoi fiecare poate fi exprimat în termeni de momente brute echivalente :

unde, în general, cum durează istoria timpului , sau numărul de puncte dacă este constantă.

Avantajul exprimării momentelor statistice în termeni de este că seturile pot fi combinate prin adăugare și nu există o limită superioară pentru valoarea .

unde îl suprascrie reprezintă istoria temporală concatenate sau combinate. Aceste valori combinate ale pot fi transformate invers în momente de origine zero (momente brute) reprezentând istoria completă a timpului concatenat

relații cunoscute între momentele de origine zero (momente brute) ( ) și momentul central ( ) sunt apoi utilizate pentru a calcula momentele centrale ale istoriei temporale legate. În cele din urmă, momentele statistice ale istoriei legate sunt calculate din momentele centrale:

Covarianță

Algoritmi foarte similari pot fi folosiți pentru a calcula covarianța . Algoritmul „banal” este:

Pentru algoritmul scris mai sus, se poate utiliza următorul cod Python:

 def național_covarianță ( date1 , date2 ):
    n = len ( data1 )
    suma12 = 0
    suma1 = suma ( date1 )
    sum2 = sum ( data2 )

    pentru i în intervalul ( n ):
        sum12 + = data1 [ i ] * data2 [ i ]

    covarianță = ( sum12 - sum1 * sum2 / n ) / n
    întoarce covarianța

În ceea ce privește variația, covarianța a două variabile aleatorii este, de asemenea, invariantă de traducere, apoi, date două valori constante Și orice, puteți scrie:

și, din nou, alegerea unei valori în intervalul de valori va stabiliza formula împotriva anulării numerice pentru ao face mai robustă împotriva sumelor mari. Luând prima valoare din fiecare set de date, algoritmul poate fi scris după cum urmează:

 def shifted_data_covariance ( dataX , dataY ):
   n = len ( dataX )
   dacă ( n < 2 ):
     returnează 0
   Kx = dataX [ 0 ]
   Ky = dataY [ 0 ]
   Ex = 0
   Ey = 0
   Exy = 0
   pentru i în intervalul ( n ):
      Ex + = dataX [ i ] - Kx
      Ey + = dataY [ i ] - Ky
      Exy + = ( dataX [ i ] - Kx ) * ( dataY [ i ] - Ky )
   return ( Exy - Ex * Ey / n ) / n

Algoritmul în doi pași calculează mai întâi eșantionul și apoi covarianța:

Algoritmul în doi pași poate fi scris după cum urmează:

 def two_pass_covariance ( data1 , data2 ):
    n = len ( data1 )

    medie1 = sumă ( date1 ) / n
    medie2 = sumă ( date2 ) / n

    covarianță = 0

    pentru i în intervalul ( n ):
        a = data1 [ i ] - medie1
        b = date2 [ i ] - medie2
        covarianță + = a * b / n

    întoarce covarianța

O versiune compensată puțin mai precisă rulează algoritmul „banal” complet pe reziduuri. Sumele finale Și ar trebui să fie zero, dar a doua trecere compensează eventualele mici erori.

O ușoară modificare a algoritmului inline pentru calcularea varianței produce un algoritm inline pentru covarianță:

 def online_covariance ( data1 , data2 ):
    medie1 = medie2 = 0.
    M12 = 0.
    n = len ( data1 )
    pentru i în intervalul ( n ):
        delta1 = ( date1 [ i ] - medie1 ) / ( i + 1 )
        mean1 + = delta1
        delta2 = ( date2 [ i ] - medie2 ) / ( i + 1 )
        mean2 + = delta2
        M12 + = i * delta1 * delta2 - M12 / ( i + 1 )
    returnează n / ( n - 1. ) * M12

Există un algoritm stabil într-un singur pas, similar cu cel văzut mai sus, care calculează co-momentul :

Asimetria aparentă din ultima ecuație se datorează faptului că , deci ambii termeni de actualizare sunt egali cu . O precizie și mai mare poate fi obținută prin prima medie și apoi folosind algoritmul rezidual stabil într-un singur pas.

Astfel, putem calcula covarianța ca

Allo stesso modo, vi è una formula per combinare le covarianze di due insiemi che può essere utilizzata per parallelizzare calcolo:

.

Note

  1. ^ a b Bo Einarsson, Accuracy and Reliability in Scientific Computing , SIAM, 1º agosto 2005, p. 47, ISBN 978-0-89871-584-2 . URL consultato il 17 febbraio 2013 .
  2. ^ a b TFChan, GH Golub and RJ LeVeque, "Algorithms for computing the sample variance: Analysis and recommendations", The American Statistician, 37 ( PDF ), 1983, pp. 242–247.
  3. ^ Nicholas Higham, Accuracy and Stability of Numerical Algorithms (2 ed) (Problem 1.10) , SIAM, 2002.
  4. ^ a b Tony F. Chan , Gene H. Golub e Randall J. LeVeque, Updating Formulae and a Pairwise Algorithm for Computing Sample Variances. ( PDF ), in Technical Report STAN-CS-79-773 , Department of Computer Science, Stanford University, 1979. .
  5. ^ Donald E. Knuth (1998). The Art of Computer Programming , volume 2: Seminumerical Algorithms , 3rd edn., p. 232. Boston: Addison-Wesley.
  6. ^ BP Welford (1962). "Note on a method for calculating corrected sums of squares and products" . Technometrics 4(3):419–420.
  7. ^ Chan, Tony F.; Golub, Gene H.; LeVeque, Randall J. (1983). Algorithms for Computing the Sample Variance: Analysis and Recommendations. The American Statistician 37, 242-247. https://www.jstor.org/stable/2683386
  8. ^ Ling, Robert F. (1974). Comparison of Several Algorithms for Computing Sample Means and Variances. Journal of the American Statistical Association, Vol. 69, No. 348, 859-866. DOI : 10.2307/2286154
  9. ^ http://www.johndcook.com/standard_deviation.html
  10. ^ DHD West (1979). Communications of the ACM , 22, 9, 532-535: Updating Mean and Variance Estimates: An Improved Method
  11. ^ Timothy B. Terriberry, Computing Higher-Order Moments Online , 2007. URL consultato il 29 aprile 2019 (archiviato dall' url originale l'8 aprile 2019) .
  12. ^ Philippe Pébay,Formulas for Robust, One-Pass Parallel Computation of Covariances and Arbitrary-Order Statistical Moments ( PDF ), in Technical Report SAND2008-6212 [ collegamento interrotto ] , Sandia National Laboratories, 2008.
  13. ^ a b Muenkeun Choi e Bert Sweetman, Efficient Calculation of Statistical Moments for Structural Health Monitoring ( PDF ), 2010. URL consultato il 17 giugno 2016 (archiviato dall' url originale il 3 marzo 2016) .

Voci correlate

Statistica Portale Statistica : accedi alle voci di Wikipedia che trattano di statistica