Algorytmy obliczania wariancji - Algorithms for calculating variance

Algorytmy obliczania wariancji odgrywają główną rolę w statystyce obliczeniowej . Kluczową trudnością w projektowaniu dobrych algorytmów dla tego problemu jest to, że wzory na wariancję mogą obejmować sumy kwadratów, co może prowadzić do niestabilności numerycznej, a także do przepełnienia arytmetycznego, gdy mamy do czynienia z dużymi wartościami.

Algorytm naiwny

Wzór na obliczenie wariancji całej populacji o rozmiarze N to:

Korzystanie korektę Bessela obliczyć bezstronnej oszacowania wariancji populacji ze skończonej próbki z n obserwacji, formuła jest:

Dlatego naiwny algorytm obliczania szacowanej wariancji ma postać:

  • Niech n ← 0, Suma ← 0, SumSq ← 0
  • Dla każdego punktu odniesienia x :
    • nn + 1
    • Suma ← Suma + x
    • SumSq ← SumSq + x × x
  • Var = (SumSq − (Suma × Suma) / n) / (n − 1)

Algorytm ten można łatwo dostosować do obliczania wariancji skończonej populacji: po prostu podziel przez N zamiast n  − 1 w ostatnim wierszu.

Ponieważ SumSq i (Sum×Sum)/ n mogą być bardzo podobnymi liczbami, anulowanie może prowadzić do tego, że precyzja wyniku będzie znacznie mniejsza niż wrodzona precyzja arytmetyki zmiennoprzecinkowej użytej do wykonania obliczeń. W związku z tym algorytm ten nie powinien być stosowany w praktyce i zaproponowano kilka alternatywnych, stabilnych numerycznie algorytmów. Jest to szczególnie złe, jeśli odchylenie standardowe jest małe w stosunku do średniej.

Obliczanie przesuniętych danych

Wariancja jest niezmienna w odniesieniu do zmian w parametrze lokalizacji , właściwości, która może być wykorzystana do uniknięcia katastrofalnego anulowania w tej formule.

z dowolną stałą, która prowadzi do nowej formuły

im bliżej wartości średniej, tym dokładniejszy będzie wynik, ale samo wybranie wartości z zakresu próbek zagwarantuje pożądaną stabilność. Jeśli wartości są małe, to nie ma problemów z sumą jego kwadratów, wręcz przeciwnie, jeśli są duże, oznacza to, że wariancja również jest duża. W każdym razie drugi termin we wzorze jest zawsze mniejszy niż pierwszy, dlatego nie może nastąpić anulowanie.

Jeśli tylko pierwsza próbka zostanie pobrana, ponieważ algorytm można napisać 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

Ta formuła ułatwia również obliczenia przyrostowe, które można wyrazić jako

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)

Algorytm dwuprzebiegowy

Alternatywne podejście, wykorzystujące inny wzór na wariancję, najpierw oblicza średnią próbki,

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

gdzie s jest odchyleniem standardowym. Podaje to następujący kod:

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 stabilny numerycznie, jeśli n jest małe. Jednak wyniki obu tych prostych algorytmów („naiwny” i „dwuprzebiegowy”) mogą nadmiernie zależeć od kolejności danych i mogą dawać słabe wyniki dla bardzo dużych zestawów danych z powodu powtarzającego się błędu zaokrąglenia w akumulacji danych. sumy. Techniki, takie jak kompensowane sumowanie, można do pewnego stopnia zwalczyć z tym błędem.

Algorytm online Welforda

Często przydatna jest możliwość obliczenia wariancji w jednym przebiegu , sprawdzając 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 dominują nad kosztami obliczeń. Dla takiego algorytm online , wykorzystując zależność nawrót jest wymagane między ilościach od której wymagane statystyki mogą być obliczone w numerycznie stabilny sposób.

Poniższe wzory mogą być użyte do aktualizacji średniej i (szacowanej) wariancji sekwencji dla dodatkowego elementu x n . Tutaj oznacza średnią próbki pierwszych n próbek , ich obciążoną wariancję próbki i ich nieobciążoną wariancję próbki .

Wzory te cierpią z powodu niestabilności liczbowej, ponieważ wielokrotnie odejmują małą liczbę od dużej liczby, która skaluje się z n . Lepszą wielkością do aktualizacji jest suma kwadratów różnic od aktualnej średniej , tutaj oznaczona :

Algorytm ten został znaleziony przez Welforda i został dokładnie przeanalizowany. Powszechne jest również oznaczanie i .

Przykładowa implementacja w Pythonie algorytmu Welforda jest podana 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
    if count < 2:
        return float("nan")
    else:
        (mean, variance, sampleVariance) = (mean, M2 / count, M2 / (count - 1))
        return (mean, variance, sampleVariance)

Ten algorytm jest znacznie mniej podatny na utratę precyzji z powodu katastrofalnego anulowania , ale może nie być tak wydajny z powodu operacji dzielenia wewnątrz pętli. W przypadku szczególnie niezawodnego dwuprzebiegowego algorytmu obliczania wariancji można najpierw obliczyć i odjąć oszacowanie średniej, a następnie użyć tego algorytmu na resztach.

Poniższy algorytm równoległy ilustruje sposób łączenia wielu zestawów statystyk obliczonych online.

Ważony algorytm przyrostowy

Algorytm można rozszerzyć o obsługę nierównych wag próbek, zastępując prosty licznik n sumą dotychczasowych wag. West (1979) sugeruje ten algorytm przyrostowy :

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)

Algorytm równoległy

Chan i in. zauważ, że opisany powyżej algorytm online Welforda jest szczególnym przypadkiem algorytmu, który działa w celu łączenia dowolnych zbiorów i :

.

Może to być przydatne, gdy na przykład wiele jednostek przetwarzających może być przypisanych do dyskretnych części danych wejściowych.

Metoda Chana do szacowania średniej jest liczbowo niestabilna, gdy i obie są duże, ponieważ błąd liczbowy w nie jest zmniejszony w taki sposób, jak w przypadku. W takich przypadkach preferuj .

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

Można to uogólnić, aby umożliwić zrównoleglenie z AVX , z procesorami graficznymi i klastrami komputerowymi oraz na kowariancję.

Przykład

Załóżmy, że wszystkie operacje zmiennoprzecinkowe używają standardowej arytmetyki podwójnej precyzji IEEE 754 . Rozważ próbkę (4, 7, 13, 16) z nieskończonej populacji. Na podstawie tej próbki szacowana średnia populacji wynosi 10, a nieobciążone oszacowanie wariancji populacji wynosi 30. Zarówno algorytm naiwny, jak i algorytm dwuprzebiegowy obliczają te wartości poprawnie.

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 dwuprzebiegowy poprawnie oblicza to oszacowanie wariancji, ale algorytm naiwny zwraca 29.333333333333332 zamiast 30.

Chociaż ta utrata precyzji może być tolerowana i postrzegana jako niewielka wada naiwnego algorytmu, dalsze zwiększanie przesunięcia powoduje, że błąd jest katastrofalny. Rozważmy próbkę ( 10 9  + 4 , 10 9  + 7 , 10 9  + 13 , 10 9  + 16 ). Ponownie szacowana wariancja populacji wynosząca 30 jest obliczana poprawnie przez algorytm dwuprzebiegowy, ale algorytm naiwny oblicza ją teraz jako -170.66666666666666. Jest to poważny problem z algorytmem naiwnym i wynika z katastrofalnego anulowania odejmowania dwóch podobnych liczb na końcowym etapie algorytmu.

Statystyki wyższego rzędu

Terriberry rozszerza formuły Chana o obliczanie trzeciego i czwartego momentu centralnego , potrzebnego na przykład przy szacowaniu skośności i kurtozy :

Tutaj znowu są sumy mocy różnic od średniej , dające

W przypadku przyrostowym (tj. ), upraszcza się to do:

Zachowując wartość , potrzebna jest tylko jedna operacja dzielenia, dzięki czemu można obliczyć statystyki wyższego rzędu przy niewielkich kosztach przyrostowych.

Przykładem algorytmu online dla kurtozy zaimplementowanego w opisany sposób jest:

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ÿ dalej rozszerza te wyniki na momenty centralne dowolnego rzędu , dla przypadków przyrostowych i parami, a następnie Pébaÿ i in. dla momentów ważonych i złożonych. Można tam również znaleźć podobne formuły na kowariancję .

Choi i Sweetman oferują dwie alternatywne metody obliczania skośności i kurtozy, z których każda może zaoszczędzić znaczne wymagania dotyczące pamięci komputera i czasu procesora w niektórych aplikacjach. Pierwsze podejście polega na obliczeniu momentów statystycznych poprzez rozdzielenie danych na przedziały, a następnie obliczenie momentów z geometrii wynikowego histogramu, który w efekcie staje się algorytmem jednoprzebiegowym dla wyższych momentów. Jedną z korzyści jest to, że obliczenia momentu statystycznego mogą być przeprowadzane z dowolną 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żna skonstruować w konwencjonalny sposób: zakres potencjalnych wartości jest podzielony na przedziały, a liczba wystąpień w każdym przedziale jest liczona i wykreślana w taki sposób, że obszar każdego prostokąta jest równy części wartości próbki w tym koszu:

gdzie i reprezentuje częstotliwość i względną częstotliwość w bin i jest całkowitym obszarem histogramu. Po tej normalizacji surowe momenty i momenty centralne można obliczyć na podstawie względnego histogramu:

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

Drugie podejście Choi i Sweetmana to analityczna metodologia łączenia momentów statystycznych z poszczególnych segmentów historii czasowej, tak aby wynikowe momenty ogólne były momentami pełnej historii czasowej. Metodologia ta może być wykorzystana do równoległego obliczania momentów statystycznych z późniejszą kombinacją tych momentów lub do kombinacji momentów statystycznych liczonych w czasie sekwencyjnym.

Jeżeli znane są zbiory momentów statystycznych: dla , to każdy z nich można wyrazić w postaci równoważnych momentów surowych:

gdzie ogólnie przyjmuje się czas trwania historii czasu lub liczbę punktów, jeśli jest stała.

Zaletą wyrażania momentów statystycznych w kategoriach , jest to, że zbiory można łączyć przez dodawanie i nie ma górnej granicy wartości .

gdzie indeks dolny reprezentuje połączoną historię czasu lub połączoną . Te połączone wartości można następnie odwrotnie przekształcić w surowe momenty reprezentujące kompletną połączoną historię czasu

Znane relacje między momentami pierwotnymi ( ) i momentami centralnymi ( ) są następnie wykorzystywane do obliczania momentów centralnych połączonej historii czasu. Wreszcie momenty statystyczne połączonej historii są obliczane z momentów centralnych:

Kowariancja

Do obliczenia kowariancji można użyć bardzo podobnych algorytmów .

Algorytm naiwny

Naiwny algorytm to:

Dla powyższego algorytmu 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 średniej

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

i ponownie wybranie wartości z zakresu wartości ustabilizuje formułę przed katastrofalnym anulowaniem, a także sprawi, że będzie ona bardziej odporna na duże sumy. Biorąc pierwszą wartość każdego zestawu danych, algorytm można zapisać jako:

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

Dwuprzebiegowy

Algorytm dwuprzebiegowy najpierw oblicza średnią próbki, a następnie kowariancję:

Algorytm dwuprzebiegowy 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ładniejsza wersja skompensowana wykonuje pełny algorytm naiwny na resztach. Ostateczne sumy i powinny wynosić zero, ale drugie przejście kompensuje każdy mały błąd.

online

Istnieje stabilny algorytm jednoprzebiegowy, podobny do algorytmu online do obliczania wariancji, który oblicza ko-moment :

Widoczna asymetria w tym ostatnim równaniu wynika z faktu, że , więc oba terminy aktualizacji są równe . Jeszcze większą dokładność można osiągnąć, najpierw obliczając środki, a następnie stosując stabilny algorytm jednoprzebiegowy na resztach.

Zatem kowariancję można obliczyć jako

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żna również dokonać niewielkiej modyfikacji, aby obliczyć kowariancję ważoną:

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 na połączenie kowariancji dwóch zestawów, który można wykorzystać do zrównoleglenia obliczeń:

Ważona wersja partiami

Istnieje również wersja ważonego algorytmu online, który wykonuje aktualizacje wsadowe: oznaczmy wagi i napisz

Kowariancję można następnie obliczyć jako

Zobacz też

Bibliografia

  1. ^ B Einarsson Bo (2005). Dokładność i niezawodność w obliczeniach naukowych . SYJAM. P. 47. Numer ISBN 978-0-89871-584-2.
  2. ^ a b c Chan, Tony F .; Golub, Gene H .; LeVeque, Randall J. (1983). „Algorytmy obliczania wariancji próbki: analiza i zalecenia” (PDF) . Statystyk amerykański . 37 (3): 242–247. doi : 10.1080/00031305.1983.10483115 . JSTOR  2683386 .
  3. ^ B c Schubert Erich; Gertz, Michael (9 lipca 2018). Numerycznie stabilne równoległe obliczanie (ko-)wariancji . ACM. P. 10. doi : 10.1145/3221269.3223036 . Numer ISBN 9781450365055. S2CID  49665540 .
  4. ^ Higham, Mikołaj (2002). Dokładność i stabilność algorytmów numerycznych (2 ed) (Problem 1.10) . SYJAM.
  5. ^ Welford, BP (1962). „Uwaga dotycząca metody obliczania skorygowanych sum kwadratów i produktów”. Technometria . 4 (3): 419–420. doi : 10.2307/1266577 . JSTOR  1266577 .
  6. ^ Donald E. Knuth (1998). The Art of Computer Programming , tom 2: Algorytmy seminumeryczne , wyd . 232. Boston: Addison-Wesley.
  7. ^ Ling, Robert F. (1974). „Porównanie kilku algorytmów do obliczania średnich próbek i wariancji”. Dziennik Amerykańskiego Towarzystwa Statystycznego . 69 (348): 859-866. doi : 10.2307/2286154 . JSTOR  2286154 .
  8. ^ http://www.johndcook.com/standard_deviation.html
  9. ^ Zachód, DHD (1979). „Aktualizacja średnich i szacunków wariancji: ulepszona metoda”. Komunikaty ACM . 22 (9): 532–535. doi : 10.1145/359146.359153 . S2CID  30671293 .
  10. ^ Chan, Tony F .; Golub, Gene H .; LeVeque, Randall J. (1979), „Aktualizacja formuł i algorytmu parami do obliczania wariancji próbki”. (PDF) , Raport techniczny STAN-CS-79-773 , Wydział Informatyki, Uniwersytet Stanforda.
  11. ^ Terriberry, Timothy B. (2007), Computing Higher-Order Moments Online , zarchiwizowane z oryginału w dniu 23 kwietnia 2014 roku , pobrane 5 maja 2008
  12. ^ Pébaÿ, Philippe (2008), „Formuły na solidne, jednoprzebiegowe równoległe obliczanie kowariancji i arbitralnej kolejności momentów statystycznych” (PDF) , Raport techniczny SAND2008-6212 , Sandia National Laboratories
  13. ^ Pebaÿ, Filip; Terriberry, Tymoteusz; Kolla, Hemanth; Bennett, Janine (2016), „Numerycznie stabilne, skalowalne formuły do ​​obliczeń równoległych i online wielowymiarowych momentów centralnych wyższego rzędu z arbitralnymi wagami” , Statystyka obliczeniowa , Springer, 31 (4): 1305–1325, doi : 10.1007/s00180- 015-0637-z , S2CID  124570169
  14. ^ B Choi Myoungkeun; Sweetman, Bert (2010), "Efficient Calculation of Statistical Moments for Structural Health Monitoring", Journal of Structural Health Monitoring , 9 (1): 13-24, doi : 10.1177/1475921709341014 , S2CID  17534100

Zewnętrzne linki