Algorytmy do obliczania wariancji - Algorithms for calculating variance


Z Wikipedii, wolnej encyklopedii

Algorytmy do obliczania wariancji odgrywają ważną rolę w statystykach obliczeniowych . Główną trudność w projektowaniu dobrych algorytmów dla tego problemu jest to, że wzory dla wariancji może obejmować sumy kwadratów, co może prowadzić do niestabilności numerycznej , a także przepełnienie gdy mamy do czynienia z dużymi wartościami.

algorytm naiwny

Wzór do obliczania wariancji całej populacji wielkości N jest:

Stosując korektę Bessela obliczyć bezstronnej oszacowania wariancji populacji od skończonego próbki z n obserwacji, formuła brzmi:

Dlatego też naiwny algorytm do obliczania szacunkowej wariancji otrzymuje brzmienie:

  • Niech n ← 0, Suma ← 0, SUMSQ ← 0
  • Dla każdego punktu odniesienia X :
    • nN + 1
    • Suma ← suma + x
    • SUMSQ ← SUMSQ + x x x
  • Var = (SUMSQ - (suma x suma) / n) / (n - 1)

Algorytm ten można łatwo przystosować do obliczenia wariancji populacji skończonej: wystarczy podzielić przez N zamiast N  - 1, na ostatniej linii.

Ze względu SUMSQ i (suma x suma) / n mogą być bardzo podobne liczby, anulowanie może prowadzić do dokładności wyniku o wiele mniej niż naturalna precyzji arytmetyki zmiennoprzecinkowej wykorzystywanego do wykonywania obliczeń. Zatem algorytm ten nie może być stosowany w praktyce, a kilka różnych numerycznie stabilny algorytmy zostały zaproponowane. Jest to szczególnie źle, jeśli odchylenie standardowe jest niewielka w stosunku do średniej. Jednak algorytm można poprawić poprzez przyjęcie sposobu zakładanej średniej .

Computing przesunięty dane

Możemy użyć własności wariancji, aby uniknąć katastrofalnych odwołania w tej formule, a mianowicie wariancja jest niezmienna w stosunku do zmian w parametrze lokalizacji

z każdej stałej, co prowadzi do nowej formule

im bliżej jest do średniej wartości bardziej dokładny wynik będzie, ale po prostu wybierając wartość wewnątrz próbki wynosić będzie gwarantować pożądaną stabilność. Jeżeli wartości są małe, to nie ma problemów z sumy kwadratów jej, wręcz przeciwnie, jeśli są duże to niekoniecznie oznacza, że wariancja jest duża, jak również. W każdym przypadku, drugi składnik we wzorze jest zawsze mniejsza niż pierwsza zatem nie mogą występować anulowania.

Jeżeli weźmiemy pod uwagę tylko pierwszą próbkę jako algorytm może być napisany w języku programowania Python jako

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

Wzór ten umożliwia również przyrostową obliczenia, które mogą być wyrażone jako

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)

Algorytm dwuprzejściowa

Alternatywne podejście, z użyciem innego wzoru, różniącym najpierw oblicza średnią przykładowego

,

a następnie oblicza się sumę kwadratów różnic od średniej,

,

gdzie s jest odchyleniem standardowym. Jest to oblicza się według następującego Pseudokod:

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

Ten algorytm jest numerycznie za stabilny, jeśli n jest mała. Jednak wyniki obu tych prostych algorytmów ( „Naive” i „Two-pass”) może zależeć nadmiernie na uporządkowanie danych i może dać słabe wyniki dla bardzo dużych zbiorów danych ze względu na powtarzające się błąd zaokrąglenia w akumulacji kwoty. Techniki takie jak wyrównaną sumowania może być stosowany do zwalczania tego błędu do pewnego stopnia.

Algorytm jest Online Welford

Często jest użyteczne, aby móc obliczyć wariancję w jednym przebiegu, kontroli każdą wartość tylko raz; Na przykład, gdy dane są gromadzone bez wystarczającej ilości pamięci, aby zachować wszystkie wartości, lub gdy koszty dostępu do pamięci dominować tych obliczeń. Dla takiego algorytmu internetowym , wykorzystując zależność nawrót jest wymagane od ilości od której wymagane statystyki mogą zostać obliczone w sposób stabilny numerycznie.

Następujące preparaty mogą być wykorzystywane do aktualizowania średnią i wariancję (szacowanej) sekwencji, dla dodatkowego elementu x n . Tutaj x n oznacza średnią próbkę z pierwszych n próbek ( x 1 , ..., x n ) y 2 n ich próbki wariancji i σ 2 n ich wariancja populacyjna.

Preparaty te cierpią na brak stabilności numerycznym, ponieważ wielokrotnie odejmować niewielką liczbę z dużej liczbie, Skale n . Lepszym ilość do aktualizowania jest sumą kwadratów różnic obecnej średniej, tu oznaczonej :

Algorytm ten został znaleziony przez Welford, i to zostało dokładnie przeanalizowane. Oczywistym jest również do określenia i .

Implementacja Pythona na przykład algorytmu Welford jest podane poniżej.

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

Algorytm ten jest znacznie mniej podatne na utratę precyzji ze względu na katastrofalną odwołania , ale może nie być tak skuteczne ze względu na podział pracy wewnątrz pętli. Dla szczególnie wytrzymałej algorytmu dwa przejścia przez obliczenie wariancji można obliczyć najpierw odjąć oszacowanie średniej, a następnie za pomocą algorytmu na pozostałości.

Algorytm równoległy poniżej przedstawia sposób łączenia wielu zestawów danych statystycznych obliczonych w Internecie.

Ważony algorytm przyrostowy

Algorytm ten może być rozszerzony na uchwyt nierównej grubości próbki, zastępując prosty licznik N , przy czym suma ciężarów widziana pory. Zachód (1979) sugeruje to przyrostowego algorytmu :

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)

algorytm równoległy

Chan i in. zwrócić uwagę, że powyższy algorytm „On-line” jest szczególnym przypadkiem od algorytmu, który pracuje dla każdej partycji próbki w zestawach , :

,

Może to być przydatne, gdy na przykład, wielokrotne Jednostki przetwarzania mogą być przypisane do oddzielnych częściach wkład.

Metoda chan szacowania średniej liczbowo jest niestabilny podczas i oba są duże, ponieważ błąd w liczbowa nie jest skalowany w dół w taki sposób, że jest w tym przypadku. W takich przypadkach wolą .

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)

Przykład

Zakładamy, że wszystkie operacje na zmiennych użyć standardowego IEEE 754 o podwójnej precyzji arytmetyki. Rozważmy próbkę (4, 7, 13, 16) z nieskończonej populacji. Na podstawie tej próby, populacja szacowana na myśli to 10, a obiektywne oszacowanie wariancji populacji jest algorytm 30. Zarówno algorytm naiwny i dwuprzejściowa obliczyć te wartości prawidłowej.

Następne pod próbkę ( 10 8  + 4 , 10 8  + 7 , 10 8  + 13 , 10 8  + 16 ), co powoduje powstanie w tym samym przybliżony wariancji, jak w pierwszej próbce. Algorytm dwa przejścia oblicza to oszacowanie wariancji prawidłowo, a algorytm powraca 29,333333333333332 naiwnych zamiast 30.

Choć utrata dokładności może być tolerowane i postrzegane jako drobne skazy algorytmu naiwnego, zwiększając przesunięcie powoduje błąd katastrofalne. Rozważmy próbkę ( 10 9  + 4 , 10 9  + 7 , 10 9  + 13 , 10 9  + 16 ). Ponownie szacowana wariancja populacji 30 jest obliczana przez algorytm poprawnie dwa przejścia, ale algorytm naiwny teraz oblicza go jako -170.66666666666666. Jest to poważny problem z algorytmu naiwnego i to ze względu na katastrofalną odwołania w odjęciem dwóch podobnych numerów w końcowej fazie algorytmu.

statystyki wyższego rzędu

Terriberry rozciąga formuły chan do obliczania trzeci i czwarty moment centralny potrzebne na przykład przy szacowaniu skośność i kurtozę :

Tu znowu są sumy uprawnień różnic od średniej , co daje

skośność:
kurtoza:

Dla przypadku przyrostowego (tj ), to upraszcza się do:

Zachowując wartość , tylko jedna operacja podział jest potrzebne i statystyki wyższego rzędu można zatem obliczyć dla małego bieżących kosztów.

Przykładem algorytmu online dla kurtozy jest realizowany zgodnie z opisem:

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

Pébaÿ rozszerza te wyniki arbitralnej rzędu centralnych momentów , na przyrostowe a parami przypadkach, a następnie Pébaÿ et al. dla ważonych i wieloskładnikowych chwil. Można też znaleźć tam podobne formuły kowariancji ,

Choi i Sweetman oferuje dwie alternatywne metody do obliczania skośność i kurtozę, z których każdy może zaoszczędzić znaczne zapotrzebowanie na pamięć i czas procesora komputera w pewnych zastosowaniach. Pierwsze podejście jest obliczenie momentów statystycznych poprzez rozdzielenie danych na pojemnikach, a następnie obliczanie momentom geometrii otrzymanego wykresu, który staje się wtedy algorytm jednego przebiegu dla wyższych momentach. Jedną z zalet jest to, że obliczenia chwila statystyczne mogą być prowadzone do arbitralnej dokładnością tak, że obliczenia mogą być dostrojone do precyzji, np formatu przechowywania danych lub oryginalnego sprzętu pomiarowego. Względny histogram zmiennej losowej może być skonstruowane w sposób konwencjonalny: zakres potencjalnych wartości dzieli się do pojemników i liczba powtórzeń w każdym przedziale są zliczane i naniesione tak, że obszar każdego prostokąta jest równa części wartości próbnych w tym bin:

gdzie i reprezentuje częstotliwość i względną częstotliwość w kosza i to łączna powierzchnia histogramu. Po normalizacji, że surowe momentów i momentów z można obliczyć względną histogramie:

gdzie górny wskazuje momenty są obliczane na podstawie histogramu. Dla stałej szerokości bin te dwa wyrażenia można uprościć za pomocą :

Drugie podejście z Choi i Sweetman jest metoda analityczna połączyć momenty statystyczne z poszczególnych segmentów czasowych historii tak, że uzyskany ogólny momenty są te całkowitego czasu historii. Ta metoda może być stosowana do równoległego obliczeń statystycznych momentów z kolejnej kombinacji tych momentach lub w połączeniu z momentów statystycznych obliczane w czasie kolejnych.

Jeśli zestawy momentów statystycznych są znane: dla , to każdy może być wyrażone w równoważnych surowych momentach:

gdzie zazwyczaj brane są za okres od czasu historii, czy liczba punktów, jeśli jest stała.

Korzyść wyrażania momentów statystycznych w zakresie jest to, że zbiory mogą być łączone przez dodanie, i nie ma górnej granicy wartości .

gdzie indeks dolny oznacza łączone czasowego historii lub w połączeniu . Połączenie tych wartości można następnie odwrotnie przekształcone surowego momentów przedstawiają kompletny łączone czasowego historię

Znane relacje między surowych momentach ( ) i centralnych momentach ( ) są następnie wykorzystywane do obliczania centralne momenty połączonego czasu historii. Wreszcie, momenty statystyczne połączonego historii są obliczane z centralnych momentach:

kowariancji

Bardzo podobne algorytmy mogą być wykorzystane do obliczenia kowariancji .

algorytm naiwny

Algorytm naiwny jest:

Dla algorytmu powyżej, można użyć następującego kodu Pythona:

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

Z oszacowaniem wartości średniej

Jak dla wariancji, kowariancja dwóch zmiennych losowych jest również przesunięcie-niezmienna, więc biorąc pod uwagę dowolne dwie wartości stałe i mogą być napisane:

i znowu wyborze wartości wewnątrz zakresu wartości ustabilizuje formułę przeciwko katastrofalnej odwołania, a także uczynić go bardziej odporne na dużych sum. Biorąc pierwszy wartość każdego zestawu danych, algorytm może być zapisany jako:

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

Dwuprzejściowa

Algorytm dwa przejścia najpierw oblicza przykładowych środków, a następnie kowariancji:

Algorytm dwuprzejściowa można zapisać jako:

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

Nieco dokładniejsze wersja skompensowana wykonuje pełen naiwnego algorytmu na reszt. Ostateczne kwoty i powinno być zero, ale drugie przejście kompensuje każdego małego błędu.

online

Stabilna algorytm jednym przebiegu istnieje, podobnie jak algorytm online dla obliczania wariancji, który oblicza współpracy chwilę :

Widoczna asymetria tego ostatniego równania wynika z faktu, że , więc oba terminy aktualizacji są równe . Jeszcze większą dokładność można osiągnąć przez pierwszy obliczanie środka, a następnie za pomocą stabilnego przejścia algorytmu na pozostałości.

W ten sposób możemy obliczyć jako kowariancji

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)

Możemy również zrobić małą modyfikację do obliczenia kowariancji ważonej:

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)

Podobnie, istnieje wzór łączenia kowariancji dwóch zestawów, które można wykorzystać do parallelize obliczenia:

Ważona wersja batched

Wersja ważonej algorytm online, które nie batched zaktualizowanej również istnieje: niech oznaczają ciężary, możemy napisać

Kowariancja może być obliczana jako

Zobacz też

Referencje

  1. ^ B Bo Einarsson (1 sierpnia 2005). Dokładność i niezawodność w Scientific Computing . Siam. p. 47. ISBN  978-0-89871-584-2 . Źródło 17 luty 2013 .
  2. ^ B T.F.Chan GH Golub RJ leveque (1983). Algorytmy obliczania wariancji próbki: Analiza i rekomendacje”, amerykański statystyk, 37" (PDF) : 242-247.
  3. ^ B Schubert Erich; Gertz, Michael (9.07.2018). „Numerycznie stabilny równolegle obliczenie wariancji (współ)” . ACM: 10. doi : 10,1145 / +3221269,3223036 . ISBN  9781450365055 .
  4. ^ Higham Nicholas (2002). Dokładność i Stabilność numeryczne algorytmów (2 ed) (problem 1,10) . Siam.
  5. ^ BP Welford (1962). „Uwaga dotycząca sposobu obliczania skorygowanych sumy kwadratów i produktów” . Technometrics 4 (3): 419-420.
  6. ^ Donald E. Knuth (1998). Sztuka programowania , tom 2: Seminumerical algorytmów ., 3. EDN, s. 232. Boston: Addison-Wesley.
  7. ^ Chan Tony F .; Golub Gene H .; Leveque Randall J. (1983). Algorytmy do obliczania wariancji próbki: Analiza i zalecenia. Amerykański Statystyk 37, 242-247. https://www.jstor.org/stable/2683386
  8. ^ Ling, Robert F. (1974). Porównanie kilku algorytmów do obliczania przykładowych średnich i wariancji. Journal of American Statistical Association, Vol. 69, nr 348, 859-866. doi : 10,2307 / 2286154
  9. ^ http://www.johndcook.com/standard_deviation.html
  10. ^ DHD Zachód (1979). Communications of the ACM , 22, 9, 532-535: Aktualizowanie średnią i wariancję wartości szacunkowych: ulepszonej metody
  11. ^ Chan Tony F ; Golub, Gene H. ; Leveque Randall J. (1979), „Aktualizacji Formuły i przez parowanie algorytm do obliczania przykładowych wariancji”. (PDF) , Technical Report STAN-CS-79-773 , Wydział Informatyki, Uniwersytet w Stanford,
  12. ^ Terriberry Timothy B. (2007), Computing wyższego rzędu Moments Online
  13. ^ Pébaÿ Philippe (2008), "Wzory do wytrzymałe, One-pass Parallel Obliczanie kowariancji i arbitralny zamówienie momentów statystycznych" (PDF) , Technical Report SAND2008-6212 , Sandia National Laboratories
  14. ^ Pébaÿ Philippe; Terriberry, Timothy; Kolla, Hemanth; Bennett, Janine (2016), „numerycznie stabilne, skalowalne i Wzory do Parallel Online Obliczenie wyższego rzędu Moments wieloczynnikowa Środkowej z dowolną Wag” , Computational Statystyki , Springer
  15. ^ B Choi Myoungkeun; Sweetman, Bert (2010), "Efektywne Obliczanie momentów statystycznych Structural Health Monitoring dla" , Journal of Structural Health Monitoring , 9 (1)

Linki zewnętrzne