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 :
- n ← n + 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
- ^ B Einarsson Bo (2005). Dokładność i niezawodność w obliczeniach naukowych . SYJAM. P. 47. Numer ISBN 978-0-89871-584-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 .
- ^ 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 .
- ^ Higham, Mikołaj (2002). Dokładność i stabilność algorytmów numerycznych (2 ed) (Problem 1.10) . SYJAM.
- ^ 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 .
- ^ Donald E. Knuth (1998). The Art of Computer Programming , tom 2: Algorytmy seminumeryczne , wyd . 232. Boston: Addison-Wesley.
- ^ 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 .
- ^ http://www.johndcook.com/standard_deviation.html
- ^ 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 .
- ^ 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.
- ^ Terriberry, Timothy B. (2007), Computing Higher-Order Moments Online , zarchiwizowane z oryginału w dniu 23 kwietnia 2014 roku , pobrane 5 maja 2008
- ^ 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
- ^ 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
- ^ 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