Algoritmy pro výpočet rozptylu - Algorithms for calculating variance


z Wikipedie, otevřené encyklopedie

Algoritmy pro výpočet rozptylu hrají významnou roli ve výpočetní statistiky . Klíčovým obtíže při navrhování dobrých algoritmů tohoto problému je to, že vzorce pro rozptyl může zahrnovat sumy čtverců, což může vést k numerické nestabilitě stejně jako přetečení při jednání s velkými hodnotami.

naivní algoritmus

Vzorec pro výpočet rozptyl celé populace o 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:

Proto je naivní algoritmus pro výpočet odhadované rozptyl je dán následující:

  • Nechť n ← 0, Sum ← 0, SUMSQ ← 0
  • Pro každou nulových bodů x :
    • nn + 1
    • Sum ← Sum + x
    • SUMSQ ← SUMSQ + x x x
  • Var = (SUMSQ - (součet x Sum) / n) / (n - 1)

Tento algoritmus lze snadno přizpůsobit pro výpočet rozptyl konečnou populaci: jednoduše dělit N namísto n  - 1 na posledním řádku.

Vzhledem k tomu, SUMSQ a (Sum × Sum) / n mohou být velmi podobná čísla, zrušení může vést k přesnosti výsledku, jehož má být mnohem nižší, než je vlastní přesnost floating-point aritmetika použité k provedení výpočtu. Proto tento algoritmus by neměl být použit v praxi, a několik střídavé, numericky stabilní, byly navrženy algoritmy. To je obzvláště špatné, pokud je směrodatná odchylka malá vzhledem k průměru. Nicméně, algoritmus může být zlepšena tím, že přijme metodu předpokládaného průměru .

Computing posunul dat

Můžeme použít vlastnost rozptylu, aby se zabránilo katastrofální zrušení v tomto vzorci totiž rozptyl je neměnný s ohledem na změny v parametru umístění

se všechny konstantní, což vede k novému vzorce

Čím blíže je střední hodnota, tím přesnější bude výsledek, ale jen vybírá hodnotu uvnitř vzorky pohybují zaručí požadovanou stabilitu. V případě, že hodnoty jsou malé pak nejsou žádné problémy se součtem jeho náměstích, na druhou stranu, pokud jsou velké to nutně znamená, že rozptyl je velký stejně. V každém případě je druhý člen ve vzorci je vždy menší než proto se mohou vyskytovat první č zrušení.

Vezmeme-li pouze první vzorek, jak algoritmus může být napsán v jazyce Python programovacím jazyce 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 usnadňuje také přírůstkové výpočet, který může být vyjádřen 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)

Dvojfázové algoritmus

Alternativní přístup, použít jiný vzorec pro rozptyl, nejprve vypočítá průměr vzorku,

,

a pak vypočítá součet čtverců rozdílů od průměru,

,

kde s je směrodatná odchylka. To je dáno následující pseudokód:

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 číselně stabilní, pokud n je malý. Nicméně, výsledky obou těchto jednoduchých algoritmů ( „Naivní“ a „Two-pass“) může záviset nadměrně na uspořádání dat a může dát špatné výsledky pro velmi velké soubory dat v důsledku opakovaného roundoff chyby v akumulaci částky. Techniky, jako je kompenzován součet může být použit v boji proti této chyby do určité míry.

Welford je online algoritmus

Často je užitečné, aby bylo možné vypočítat rozptyl v jednom průchodu, kontrola každou hodnotu pouze jednou; Například, pokud jsou údaje shromažďovány, aniž by dostatek místa pro uložení, aby všechny hodnoty, nebo když náklady na přístup do paměti dominují ty výpočtu. 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é.

Tyto vzorce mohou být použity k aktualizaci střední a (odhad) rozptyl sekvence, o další prvek x n . Zde x n značí střední hodnoty výběru prvních n vzorků ( x 1 , ..., x n ), to 2 n jejich výběrový rozptyl, a σ 2 n jejich populace rozptyl.

Tyto vzorce trpí numerické nestability, protože opakovaně odečíst malé množství z velké číslo, které váhy s n . Lepší množství pro aktualizaci je součet čtverců rozdílů z aktuální průměr, zde označený :

Tento algoritmus byl nalezen Welfordu, a to byl důkladně analyzován. To je také obyčejné pro označení a .

Implementace příklad Python pro Welford algoritmu je uveden níže.

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

Tento algoritmus je mnohem méně náchylný ke ztrátě přesnosti díky katastrofálnímu zrušení , ale nemusí být tak efektivní, protože operace dělení uvnitř smyčky. Pro obzvláště robustní Dvojfázové algoritmus pro výpočet rozptylu, lze vypočítat první a odečíst odhad průměrné hodnoty, a pak použít tento algoritmus na rezidua.

Paralelní algoritmus níže ukazuje, jak sloučit více sad statistik vypočítaných online.

Vážený inkrementální algoritmus

Algoritmus může být rozšířen pro zpracování nestejné hmotnosti vzorků, nahrazující jednoduché čítače n se součtem hmotností dosud viděl. West (1979) naznačuje, tato inkrementální algoritmus :

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)

paralelní algoritmus

Chan a kol. Všimněte si, že výše uvedené „On-line“ algoritmus je zvláštní případ algoritmu, který pracuje pro jakékoliv rozdělení vzorku do sad , :

,

To může být užitečné, když, například, mohou být vícenásobné jednotky pro zpracování přiřazen diskrétních částí vstupu.

Chan metoda odhadu průměr je numericky nestabilní, kdy a jak jsou velké, protože číselná chyba není zmenšen tak, že je to v případě. V takových případech přednost .

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)

Příklad

Předpokládají, že všechny operace plovoucí čárkou používat standardní IEEE 754 dvojitou přesností aritmetiky. Předpokládejme vzorek (4, 7, 13, 16) z nekonečnou populaci. Na tomto vzorku založené odhadovaná populace průměr je 10, a nestranný odhad populačního rozptylu je 30. Jak naivní algoritmus a Dvojfázové Algoritmus výpočtu těchto hodnot správně.

Dále v úvahu vzorku ( 10 8  + 4 , 10 8  + 7 , 10 8  + 13 , 10 8  + 16 ), který vede ke stejným odhadovaný rozptyl jako první vzorek. Dvojfázové algoritmus počítá tento rozptyl odhadu správně, ale naivní algoritmus vrací 29,333333333333332 namísto 30.

I když tato ztráta přesnosti může být snesitelný a zobrazil jako menší chybu z naivního algoritmu, což dále zvyšuje posun dělá chyba katastrofální. Předpokládejme vzorek ( 10 9  + 4 , 10 9  + 7 , 10 9  + 13 , 10 9  + 16 ). Opět odhadnutá populace rozptyl 30 je vypočítána správně Dvojfázové algoritmus, ale naivní algoritmus nyní počítá jako -170.66666666666666. Jedná se o závažný problém s naivní algoritmus a je kvůli katastrofálním zrušení v odečtením dvou podobných čísel v konečné fázi algoritmu.

Statistiky vyššího řádu

Terriberry rozšiřuje vzorců Chan, aby výpočet třetí a čtvrté centrální momenty , potřebné například při odhadování šikmost a špičatost :

Zde jsou opět sumy sil rozdíly od průměru , přičemž

šikmost:
špičatost:

V případě inkrementálního případě (tj ), to zjednoduší na:

Zachováním hodnotu , je potřeba pouze jedna operace dělení a statistika vyššího řádu tak může být vypočítána za málo přírůstkových nákladů.

Příkladem online algoritmus pro špičatost realizován, jak bylo 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

    kurtosis = (n*M4) / (M2*M2) - 3
    return kurtosis

Pébaÿ dále rozšiřuje tyto výsledky na libovolný řádu centrální momenty , pro přírůstkové a pairwise případech, a následně Pébaÿ et al. pro vážené a složených momentů. Dá se tam také najít podobné vzorce pro kovariance .

Choi a Sweetman nabízejí dvě alternativní metody pro výpočet šikmost a špičatost, z nichž každý může ukládat požadavky podstatné počítačové paměti a procesorového času v některých aplikacích. První přístup je výpočet statistické momenty oddělením data do koše a vypočítávání momenty z geometrie výsledného histogramu, který účinně stane algoritmus jednofázové pro vyšší momenty. Jednou z výhod je, že statistické výpočty moment se může provádět na libovolném přesností tak, že výpočty mohou být naladěn na přesnost, například, formát pro ukládání dat, nebo původní měřicí hardware. Relativní Histogram náhodné proměnné mohou být vytvořeny běžným způsobem: rozsah možných hodnot, je rozdělena do zásobníků a počet výskytů v každém zásobníku se počítají a vynesou tak, že plocha každé obdélníku se rovná část hodnot vzorků v rámci tohoto koše:

kde a představuje frekvence a relativní četnost v zásobníku a je celková plocha histogramu. Po této normalizace se surové momenty a centrální momenty lze vypočítat z relativního histogramu:

kde index označuje momenty jsou vypočteny z histogramu. Pro konstantní šířce přihrádky mohou být tyto dva výrazy zjednodušené použití :

Druhý přístup z Choi a Sweetman je analytická metoda kombinuje statistické momenty z jednotlivých segmentů časového historii tak, že výsledné celkové momenty jsou úplné časové historie. Tento postup by mohl být použit pro paralelního výpočtu statistických momentů s následnou kombinací těchto momentů, nebo kombinaci statistických momentů vypočtených v po sobě jdoucích krát.

Pokud jsou známy sady statistických momentů: pro , pak každý může být vyjádřena z hlediska ekvivalentní surových momentů:

kde je obecně vzat být délka časového historii, nebo počet bodů, pokud je konstantní.

Přínosem vyjadřující statistické momenty z hlediska je, že sady mohou být kombinovány přidáním, a neexistuje žádný horní limit na hodnotě .

kde index reprezentuje zřetězené časové historii nebo kombinované . Tyto kombinované hodnoty pak mohou být nepřímo přeměněn surových momenty, které představují kompletní zřetězené časové historii

Známé vztahy mezi surových momentů ( ) a centrální momenty ( ) jsou pak použity k výpočtu centrální momenty zřetězené časové historie. A konečně, statistické momenty zřetězené historie jsou počítány z centrálních momentů:

kovariance

Velmi podobné algoritmy mohou být použity pro výpočet kovariance .

naivní algoritmus

Naivní algoritmus je:

U algoritmu výše, by se dalo použít následující kód v jazyce Python:

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 rozptylu, kovariance dvou náhodných proměnných je také posun-neměnný, takže vzhledem k tomu, žádné dvě konstantní hodnoty , a to může být psáno:

a znovu výběru hodnoty uvnitř rozsahu hodnot bude stabilizovat formuli proti katastrofickým zrušení stejně jako aby bylo více odolné vůči velkých částek. Vezmeme-li první hodnotu každého souboru dat, algoritmus lze zapsat 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

Dvojfázové

Dvojfázové algoritmus nejprve vypočítá ukázkové prostředky, a pak kovariance:

Dvojfázové 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

O něco přesnější kompenzován verze provádí plnou naivní algoritmus na rezidua. Konečné částky a měla být nulová, ale ve druhém průchodu kompenzuje malé chyby.

Online

Stabilní jednoprůchodový algoritmus existuje, podobně jako on-line algoritmus pro výpočet rozptylu, který počítá spolupráce moment :

Zdánlivá asymetrie v tom posledním rovnice je vzhledem k tomu, že , takže obě aktualizace termíny jsou rovné . Ještě větší přesnost může být dosaženo tím, že nejprve výpočetní prostředky, pak pomocí stabilní jednofázové algoritmus na zbytky.

Můžeme tedy počítat kovariance 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)

Můžeme také udělat malou změnu k výpočtu váženého 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)

Stejně tak je zde vzorec pro kombinování kovariance dvou sad, které mohou být použity k paralelizovat výpočet:

Vážený dávkovaný verze

Verze váženého on-line algoritmem, který se porcovaný také aktualizovaný existuje nechť označují závaží, můžeme napsat

Kovariance může být potom vypočítán jako

viz též

Reference

  1. ^ B Bo Einarsson (01.8.2005). Přesnost a spolehlivost v Scientific Computing . SIAM. p. 47. ISBN  978-0-89871-584-2 . Vyvolány 17 February rok 2013 .
  2. ^ B T.F.Chan, GH Golub a RJ LEVEQUE (1983). " " Algoritmy pro výpočet výběrový rozptyl: Analýza a doporučení“, The American statistik, 37" (PDF) : 242-247.
  3. ^ B Schubert, Erich; Gertz, Michael (07.9.2018). „Číslicově stabilní paralelní výpočet (spolu) rozptylu“ . ACM: 10. doi : 10,1145 / 3.221.269,3223036 . ISBN  9781450365055 .
  4. ^ Higham, Nicholas (2002). Přesnost a stabilita numerických algoritmů (2 ED) (problém 1.10) . SIAM.
  5. ^ BP Welford (1962). „Poznámka k metodě pro výpočet korigovaných sumy čtverců a produkty“ . Technometrics 4 (3): 419-420.
  6. ^ Donald E. Knuth (1998). Umění programování , Volume 2: Seminumerical algoritmech ., 3. EDN, str. 232. Boston: Addison-Wesley.
  7. ^ Chan, Tony F .; Golub, Gene H .; LEVEQUE, Randall, J. (1983). Algoritmy pro výpočet výběrový rozptyl: analýza a doporučení. Americký Statistician 37, 242-247. https://www.jstor.org/stable/2683386
  8. ^ Ling, Robert F. (1974). Porovnání několika Algoritmy pro výpočet vzorkem se rozumí a rozptyl. Journal of American statistické asociace, Vol. 69, No. 348, 859-866. doi : 10,2307 / 2286154
  9. ^ http://www.johndcook.com/standard_deviation.html
  10. ^ DHD West (1979). Komunikace ACM , 22, 9, 532-535: Aktualizace střední hodnoty a rozptylu odhadů: vylepšený způsob
  11. ^ Chan, Tony F. ; Golub, Gene H. ; LEVEQUE, Randall J. (1979), "Aktualizace Vzorce a Pairwise algoritmus pro výpočet ukázkové odchylky." (PDF) , Technical Report STAN-CS-79-773 , ministerstvo informatiky, Stanford University,
  12. ^ Terriberry, Timothy B. (2007), Computing vyššího řádu Moments Online
  13. ^ Pébaÿ, Philippe (2008), "Vzorce pro Robust, jednoprůchodový paralelního výpočtu kovariance a svévolné-order Statistical Moments" (PDF) , technická zpráva SAND2008-6212 , Sandia National Laboratories
  14. ^ Pébaÿ, Philippe; Terriberry, Timothy; Kolla, Hemanth; Bennett, Janine (2016), „číslicově stabilní, škálovatelné Vzorce pro paralelní a online výpočtu vyššího řádu Multivariate centrální momenty s libovolnou Hmotnosti“ , delšími Statistika , Springer
  15. ^ B Choi, Myoungkeun; Sweetman, Bert (2010), "Efektivní Výpočet statistických momentů pro monitorování Structural Health" , Journal of Structural Health Monitoring , 9 (1)

externí odkazy