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 :
    • nn + 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

  1. ^ 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.
  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 .
  3. ^ 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 .
  4. ^ Higham, Nicholas (2002). Přesnost a stabilita numerických algoritmů (2. vydání) (problém 1.10) . SIAM.
  5. ^ 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 .
  6. ^ Donald E. Knuth (1998). The Art of Computer Programming , volume 2: Seminumerical Algorithms , 3rd edn., S. 232. Boston: Addison-Wesley.
  7. ^ 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 .
  8. ^ http://www.johndcook.com/standard_deviation.html
  9. ^ 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 .
  10. ^ 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.
  11. ^ Terriberry, Timothy B. (2007), Computing Higher-Order Moments Online , archivováno z originálu dne 23. dubna 2014 , vyvoláno 5. května 2008
  12. ^ 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
  13. ^ 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
  14. ^ 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

externí odkazy