Algoritmy pro výpočet rozptylu - Algorithms for calculating variance
Algoritmy pro výpočet rozptylu hrají hlavní roli ve výpočetní statistice . Klíčovým problémem při navrhování dobrých algoritmů pro tento problém je, že vzorce pro rozptyl mohou zahrnovat součty čtverců, což může vést k numerické nestabilitě a také k aritmetickému přetečení při práci s velkými hodnotami.
Naivní algoritmus
Vzorec pro výpočet rozptylu celé populace velikosti N je:
Použití Bessel je korekce pro výpočet nezkreslený odhad rozptylu populace z konečného vzorku z n pozorování, vzorec je:
Naivní algoritmus pro výpočet odhadovaného rozptylu je tedy dán následujícím:
- Nechť n ← 0, součet ← 0, SumSq ← 0
- Pro každý vztažný bod x :
- n ← n + 1
- Součet ← Součet + x
- SumSq ← SumSq + x × x
- Var = (SumSq - (Sum × Sum) / n) / (n - 1)
Tento algoritmus lze snadno upravit pro výpočet rozptylu konečné populace: jednoduše vydělte N místo n - 1 na posledním řádku.
Protože SumSq a (Sum × Sum)/ n mohou být velmi podobná čísla, zrušení může vést k tomu, že přesnost výsledku bude mnohem menší než inherentní přesnost aritmetiky s pohyblivou řádovou čárkou použitou k provedení výpočtu. Tento algoritmus by proto neměl být v praxi používán a bylo navrženo několik alternativních, numericky stabilních algoritmů. To je obzvláště špatné, pokud je standardní odchylka vzhledem k průměru malá.
Výpočet posunutých dat
Rozptyl je neměnný, pokud jde o změny v parametru umístění , což je vlastnost, kterou lze v tomto vzorci zabránit katastrofickému zrušení.
s libovolnou konstantou, což vede k novému vzorci
čím blíže je střední hodnota, tím přesnější bude výsledek, ale pouhá volba hodnoty v rozsahu vzorků zaručí požadovanou stabilitu. Pokud jsou hodnoty malé, pak nejsou problémy se součtem jeho druhých mocnin, naopak, pokud jsou velké, nutně to znamená, že i rozptyl je velký. V každém případě je druhý výraz ve vzorci vždy menší než ten první, proto nemůže dojít ke zrušení.
Pokud je odebrán pouze první vzorek, může být algoritmus zapsán v programovacím jazyce 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
Tento vzorec také usnadňuje přírůstkový výpočet, který lze vyjádřit 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)
Dvouprůchodový algoritmus
Alternativní přístup s použitím odlišného vzorce pro rozptyl nejprve vypočítá průměr vzorku,
a poté vypočítá součet druhých mocnin rozdílů od průměru,
kde s je standardní odchylka. To je dáno následujícím kódem:
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
Tento algoritmus je numericky stabilní, pokud n je malé. Výsledky obou těchto jednoduchých algoritmů („naivní“ a „dvouprůchodový“) však mohou nadměrně záviset na uspořádání dat a mohou poskytovat špatné výsledky pro velmi velké soubory dat kvůli opakované chybě zaokrouhlování při akumulaci částky. K boji s touto chybou lze do určité míry použít techniky, jako je kompenzované součty .
Welfordův online algoritmus
Často je užitečné mít možnost vypočítat rozptyl v jediném průchodu a zkontrolovat každou hodnotu pouze jednou; například když jsou data shromažďována bez dostatečného úložiště, aby byly zachovány všechny hodnoty, nebo když náklady na přístup k paměti dominují nákladům na výpočet. Pro takové on-line algoritmu , je vztah opakování mezi množství, ze kterých mohou být požadované statistické údaje vypočtené v numericky stabilním způsobem je nutné.
Následující vzorce lze použít k aktualizaci průměrné a (odhadované) odchylky posloupnosti pro další prvek x n . Zde označuje vzorek průměr prvních n vzorků , jejich předepjat vzorek rozptylu a jejich nezkreslené vzorek rozptylu .
Tyto vzorce trpí numerickou nestabilitou, protože opakovaně odečítají malé číslo od velkého čísla, které se škáluje s n . Lepší množství pro aktualizaci je součet druhých mocnin rozdílů od aktuálního průměru , zde označeného :
Tento algoritmus našel Welford a byl důkladně analyzován. Je také běžné označovat a .
Níže je uveden příklad implementace Pythonu pro Welfordův algoritmus.
# 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)
Tento algoritmus je mnohem méně náchylný ke ztrátě přesnosti v důsledku katastrofického zrušení , ale nemusí být tak účinný kvůli operaci dělení uvnitř smyčky. U obzvláště robustního dvouprůchodového algoritmu pro výpočet rozptylu lze nejprve vypočítat a odečíst odhad průměru a poté tento algoritmus použít na zbytky.
Níže uvedený paralelní algoritmus ukazuje, jak sloučit více sad statistik vypočítaných online.
Vážený inkrementální algoritmus
Algoritmus lze rozšířit tak, aby zpracovával nestejné hmotnosti vzorků, nahrazením jednoduchého čítače n součtem dosud viděných hmotností. West (1979) navrhuje tento přírůstkový algoritmus :
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)
Paralelní algoritmus
Chan a kol. Všimněte si, že výše popsaný online algoritmus společnosti Welford je speciální případ algoritmu, který funguje pro kombinování libovolných sad a :
- .
To může být užitečné, když například jednotlivým částem vstupu může být přiřazeno více procesorových jednotek.
Chanova metoda odhadu průměru je numericky nestabilní, když jsou obě velké, protože numerická chyba v není zmenšena tak, jak je tomu v případě. V takových případech raději .
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
To lze zobecnit, aby bylo možné paralelizaci s AVX , s GPU a počítačovými klastry a kovarianci.
Příklad
Předpokládejme, že všechny operace s plovoucí desetinnou čárkou používají standardní aritmetiku dvojité přesnosti IEEE 754 . Zvažte vzorek (4, 7, 13, 16) z nekonečné populace. Na základě tohoto vzorku je odhadovaný průměr populace 10 a nezaujatý odhad rozptylu populace je 30. Naivní i dvouprůchodový algoritmus tyto hodnoty vypočítávají správně.
Dále zvažte vzorek ( 10 8 + 4 , 10 8 + 7 , 10 8 + 13 , 10 8 + 16 ), který vede ke stejné odhadované variaci jako první vzorek. Dvouprůchodový algoritmus vypočítá tento odhad rozptylu správně, ale naivní algoritmus místo 29 vrátí 29,33333333333333332.
I když tato ztráta přesnosti může být tolerovatelná a považována za menší vadu naivního algoritmu, další zvýšení offsetu činí chybu katastrofickou. Zvažte vzorek ( 10 9 + 4 , 10 9 + 7 , 10 9 + 13 , 10 9 + 16 ). Odhadovaný rozptyl populace 30 je opět správně vypočítán dvouprůchodovým algoritmem, ale naivní algoritmus jej nyní vypočítá jako −170,66666666666666. Toto je vážný problém s naivním algoritmem a je důsledkem katastrofického zrušení odečtením dvou podobných čísel v konečné fázi algoritmu.
Statistiky vyššího řádu
Terriberry rozšiřuje Chanovy vzorce na výpočet třetího a čtvrtého centrálního momentu , potřebných například při odhadu šikmosti a kurtózy :
Zde jsou opět součty sil rozdílů od průměru , dávání
Pro přírůstkový případ (tj. ) To zjednodušuje:
Zachováním hodnoty je zapotřebí pouze jedna operace dělení a statistiky vyššího řádu tak lze vypočítat za malé přírůstkové náklady.
Příklad online algoritmu pro kurtózu implementovaného, jak je popsáno, je:
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ÿ dále rozšiřuje tyto výsledky na centrální momenty libovolného pořadí pro přírůstkové a párové případy a následně Pébaÿ et al. pro vážené a složené momenty. Lze zde také najít podobné vzorce pro kovarianci .
Choi a Sweetman nabízejí dvě alternativní metody pro výpočet šikmosti a špičatosti, z nichž každá může v určitých aplikacích ušetřit značné nároky na paměť počítače a čas procesoru. Prvním přístupem je vypočítat statistické momenty tak, že data rozdělíte do přihrádek a poté vypočítáte momenty z geometrie výsledného histogramu, který se ve skutečnosti stává jednoprůchodovým algoritmem pro vyšší momenty. Jednou výhodou je, že výpočty statistických momentů lze provádět s libovolnou přesností, takže výpočty lze naladit na přesnost např. Formátu pro ukládání dat nebo původního měřicího hardwaru. Relativní histogram náhodné veličiny lze sestrojit běžným způsobem: rozsah potenciálních hodnot je rozdělen do zásobníků a počet výskytů v každém zásobníku je spočítán a vykreslen tak, že plocha každého obdélníku se rovná části hodnot vzorku. v tom koši:
kde a představují frekvenci a relativní frekvenci v bin a je celkovou plochou histogramu. Po této normalizaci lze surové momenty a centrální momenty vypočítat z relativního histogramu:
kde horní index označuje, že momenty jsou vypočítány z histogramu. Pro konstantní šířku přihrádky lze tyto dva výrazy zjednodušit pomocí :
Druhý přístup od Choi a Sweetmana je analytická metodologie, která kombinuje statistické momenty z jednotlivých segmentů časové historie tak, že výsledné celkové momenty jsou momenty celé časové historie. Tuto metodiku lze použít pro paralelní výpočet statistických momentů s následnou kombinací těchto momentů nebo pro kombinaci statistických momentů počítaných v sekvenčních časech.
Pokud jsou známy sady statistických momentů: for , pak každý může být vyjádřen pomocí ekvivalentních surových momentů:
kde se obecně považuje doba trvání časové historie nebo počet bodů, je- li konstantní.
Výhodou vyjádření statistických momentů z hlediska je, že sady lze kombinovat sčítáním a neexistuje žádná horní hranice hodnoty .
kde dolní index představuje zřetězenou časovou historii nebo kombinovanou . Tyto kombinované hodnoty lze potom nepřímo transformovat do surových momentů představujících kompletní zřetězenou časovou historii
Známé vztahy mezi surovými momenty ( ) a centrálními momenty ( ) se pak použijí k výpočtu centrálních momentů zřetězené časové historie. Nakonec jsou statistické momenty zřetězené historie vypočítány z centrálních momentů:
Kovariance
K výpočtu kovariance lze použít velmi podobné algoritmy .
Naivní algoritmus
Naivní algoritmus je:
Pro výše uvedený algoritmus lze použít následující kód Pythonu:
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
S odhadem průměru
Pokud jde o rozptyl, kovariance dvou náhodných proměnných je také posunově invariantní, takže vzhledem k jakýmkoli dvěma konstantním hodnotám a lze ji zapsat:
a opět výběr hodnoty uvnitř rozsahu hodnot stabilizuje vzorec proti katastrofickému zrušení a také bude odolnější vůči velkým částkám. Vezmeme -li první hodnotu každé sady dat, algoritmus lze zapsat 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
Dvouprůchodový
Dvouprůchodový algoritmus nejprve vypočítá ukázkové prostředky a poté kovarianci:
Dvouprůchodový algoritmus může být zapsán 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
Mírně přesnější kompenzovaná verze provádí plně naivní algoritmus zbytků. Konečné součty a měly by být nulové, ale druhý průchod kompenzuje jakoukoli malou chybu.
Online
Existuje stabilní jednoprůchodový algoritmus, podobný online algoritmu pro výpočet rozptylu, který vypočítává souběžný moment :
Zjevná asymetrie v této poslední rovnici je dána skutečností, že obě podmínky aktualizace jsou stejné . Ještě větší přesnosti lze dosáhnout nejprve výpočtem prostředků a poté použitím stabilního jednoprůchodového algoritmu na zbytky.
Kovarianci lze tedy vypočítat 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)
Lze také provést malou úpravu pro výpočet vážené kovariance:
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)
Podobně existuje vzorec pro kombinaci kovariancí dvou sad, které lze použít k paralelizaci výpočtu:
Vážená dávková verze
Existuje také verze váženého online algoritmu, který provádí dávkovou aktualizaci: označte váhy a zapište
Kovarianci lze pak vypočítat jako
Viz také
Reference
- ^ a b Einarsson, Bo (2005). Přesnost a spolehlivost ve vědeckých počítačích . SIAM. p. 47. ISBN 978-0-89871-584-2.
- ^ a b c Chan, Tony F .; Golub, Gene H .; LeVeque, Randall J. (1983). „Algoritmy pro výpočet rozptylu vzorků: Analýza a doporučení“ (PDF) . Americký statistik . 37 (3): 242–247. doi : 10,1080/00031305.1983.10483115 . JSTOR 2683386 .
- ^ a b c Schubert, Erich; Gertz, Michael (9. července 2018). Numericky stabilní paralelní výpočet (ko-) rozptylu . ACM. p. 10. doi : 10,1145/3221269.3223036 . ISBN 9781450365055. S2CID 49665540 .
- ^ Higham, Nicholas (2002). Přesnost a stabilita numerických algoritmů (2. vydání) (problém 1.10) . SIAM.
- ^ Welford, BP (1962). „Poznámka k metodě výpočtu opravených součtů čtverců a produktů“. Technometrics . 4 (3): 419–420. doi : 10,2307/1266577 . JSTOR 1266577 .
- ^ Donald E. Knuth (1998). The Art of Computer Programming , volume 2: Seminumerical Algorithms , 3rd edn., S. 232. Boston: Addison-Wesley.
- ^ Ling, Robert F. (1974). „Porovnání několika algoritmů pro výpočet vzorových prostředků a variant“. Journal of the American Statistical Association . 69 (348): 859–866. doi : 10,2307/2286154 . JSTOR 2286154 .
- ^ http://www.johndcook.com/standard_deviation.html
- ^ West, DHD (1979). „Aktualizace průměrů a odhadů odchylek: vylepšená metoda“. Komunikace ACM . 22 (9): 532–535. doi : 10,1145/359146,359153 . S2CID 30671293 .
- ^ Chan, Tony F .; Golub, Gene H .; LeVeque, Randall J. (1979), „Aktualizace vzorců a párový algoritmus pro výpočet vzorových odchylek“. (PDF) , Technická zpráva STAN-CS-79-773 , Katedra informatiky, Stanfordská univerzita.
- ^ Terriberry, Timothy B. (2007), Computing Higher-Order Moments Online , archivováno z originálu dne 23. dubna 2014 , vyvoláno 5. května 2008
- ^ Pébaÿ, Philippe (2008), „Vzorce pro robustní, jednoprůchodové paralelní výpočty kovariancí a statistických momentů libovolného řádu“ (PDF) , technická zpráva SAND2008-6212 , Sandia National Laboratories
- ^ Pébaÿ, Philippe; Terriberry, Timothy; Kolla, Hemanth; Bennett, Janine (2016), „Numerically Statable, Scalable Formulas for Parallel and Online Computation of Higher-Order Multivariate Central Moments with Arbitrary Weights“ , Computational Statistics , Springer, 31 (4): 1305–1325, doi : 10.1007/s00180- 015-0637-z , S2CID 124570169
- ^ a 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