Zapište =NORM.DIST(115,100,15,TRUE) do buňky a Excel vrátí hodnotu 0.8413447 bez jakýchkoli okolků. Toto volání vypadá jako pouhé vyhledání v tabulce. Ale není tomu tak. Za tímto jediným číslem se skrývá kumulativní normální rozdělení, integrál bez uzavřené formy, a za funkcemi CHISQ.INV.RT a BETA.DIST stojí speciální funkce, které musí pečlivá knihovna skutečně vyhodnocovat a nikoli pouze odhadovat. Tabulkový komponent, který deklaruje kompatibilitu s Excelem, musí tyto hodnoty reprodukovat do posledního čísla, které Excel zobrazuje, což znamená implementovat skutečné numerické metody a ne pouze názvy funkcí
Knihovna HotXLS implementuje více než padesát těchto statistických funkcí a práce, díky níž jsou správné, je z řádku vzorců téměř neviditelná. Toto je průvodce tím, jak je výpočetní jádro počítá: sdílené jádro speciálních funkcí, rozhodování o větvení, které udržuje aritmetiku stabilní, a jedna chyba inverzního normálního rozdělení, která se dlouho skrývala v okrajové části (tail), protože běžné případy se nefunkčního řádku kódu nikdy nedotkly
Jedno volání v listu, padesát rozdělení za ním
Tyto funkce pokrývají skupiny, po kterých statistický sešit sahá. Patří sem normální rozdělení, NORM.DIST a NORM.S.DIST s jejich inverzemi; rodina rozdělení gamma a chí-kvadrát, GAMMA.DIST, CHISQ.DIST, CHISQ.DIST.RT, CHISQ.INV.RT; rozdělení beta, BETA.DIST a BETA.INV; výběrová rozdělení T.DIST, T.DIST.2T, F.DIST a F.INV; diskrétní dvojice BINOM.DIST a POISSON.DIST; a pomocné statistické funkce jako CONFIDENCE.T a CONFIDENCE.NORM. Z pohledu volajícího je každá z nich jediným vzorcem. Nastavíte vstupy v buňkách, necháte sešit provést výpočet a přečtete výsledek
var
wb: IXLSWorkbook;
sh: IXLSWorksheet;
begin
wb := TXLSWorkbook.Create;
sh := wb.Sheets.Add;
sh.Range['A1', 'A1'].Value := 115; // observation
sh.Range['A2', 'A2'].Value := 100; // mean
sh.Range['A3', 'A3'].Value := 15; // standard deviation
// The XLS formula parser uses ';' as the argument separator.
Writeln(wb.Calculate('=NORM.DIST(A1;A2;A3;TRUE())')); // 0.8413447
Writeln(wb.Calculate('=CHISQ.INV.RT(0.05;10)')); // 18.3070381
Writeln(wb.Calculate('=BETA.DIST(0.5;2;3;TRUE())')); // 0.6875
end;Metoda Calculate v sešitu kompiluje a vyhodnocuje ad-hoc vzorec vůči aktivnímu listu a vrací typ Variant. Jeden detail může napoprvé uživatele zmást: parser vzorců na pozadí Calculate přijímá jako oddělovač argumentů středník, takže se píše =SUM(A1;B1) a nikoli =SUM(A1,B1). Vzorce uložené v buňkách si však ponechávají čárku standardní pro Excel. Stejný vyhodnocovač odbavuje každou z níže popsaných statistických funkcí, takže jakmile jedna z nich funguje v Calculate, ostatní následují stejnou cestu
Dvě funkce, na kterých staví vše ostatní
Většina kumulativních rozdělení v této sadě se nepočítá sčítáním nebo integrováním jejich vlastních definic. Počítají se ze dvou speciálních funkcí: regularizované dolní neúplné funkce gamma, zapisované jako P(a, x), a regularizované neúplné funkce beta, zapisované jako Ix(a, b). Interně se jedná o pomocné funkce, o které se opírají příslušné vyhodnocovače, a tento řetězec je krátký. Kumulativní distribuční funkce (CDF) chí-kvadrát je CDF rozdělení gamma s parametrem tvaru df/2 a měřítkem 2. CDF rozdělení gamma je přímo P(a, x). Kumulativní funkce t, F a binomického rozdělení jsou hodnoty regularizované neúplné funkce beta při odpovídajících argumentech. CDF Poissonova rozdělení je horní neúplná funkce gamma Q. Implementujte funkce gamma a beta správně a tucet rozdělení zdědí jejich přesnost zdarma
Slovo „regularizovaná“ (regularized) je zde klíčové. Původní neúplná funkce gamma roste jako faktoriál a původní integrál funkce beta může způsobit podtečení (underflow) nebo přetečení (overflow) dlouho předtím, než se dospěje k výsledku. Regularizované formy jsou poděleny úplnou funkcí gamma nebo beta, takže se pohybují výhradně v intervalu od nuly do jedné, což je přesně rozsah, který zaujímá pravděpodobnost. Tato normalizace umožňuje, aby stejná procedura sloužila pro rozdělení chí-kvadrát se dvěma stupni volnosti i se dvěma sty stupni volnosti, aniž by mezilehlé členy přetekly limit typu double. Vysvětluje to také, proč nepočítáte CDF sčítáním dlouhé řady členů hustoty: každý člen nese svou vlastní zaokrouhlovací chybu, chyby se v průběhu řady sčítají a regularizovaná speciální funkce se této sumě zcela vyhýbá tím, že namísto toho vyhodnocuje rychle konvergující řadu nebo řetězový zlomek
Řada pod diagonálou, řetězový zlomek nad ní
Procedura neúplné funkce gamma činí před výpočtem jedno rozhodnutí: porovnává hodnotu x s hodnotou a + 1. Tato hranice není náhodná. Rozvoj mocninné řady funkce P(a, x) konverguje rychle, když je x malé vůči a, a pomalu, případně nepoužitelně, když je x velké. Řetězový zlomek má opačnou vlastnost. Jádro proto používá mocninnou řadu pro x menší než a + 1 a Lentzův řetězový zlomek pro x rovno nebo větší než a + 1, přičemž každá větev provádí pouze tu práci, pro kterou je vhodná
Řetězový zlomek vyžaduje jedno zabezpečení. Lentzova metoda pracuje s průběžným čitatelem a jmenovatelem a invertuje jmenovatel v každém kroku. Pokud se některý z nich přiblíží nule, inverze selže. Nápravou je minimální limit: kdykoli mezilehlý člen klesne pod hodnotu zhruba 1e-30, je na tuto hodnotu oříznut. To udržuje rekurenci konečnou, aniž by to narušilo konvergovanou hodnotu. Stejné oříznutí se ze stejného důvodu objevuje i v řetězovém zlomku neúplné funkce beta. Jde o malou konstantu plnící důležitou roli — představuje rozdíl mezi stabilním vyhodnocením a dělením hodnotou nerozeznatelnou od nuly
Horní okrajová část, Q(a, x), je jednoduše 1 minus P(a, x), a takto se počítá kumulativní větev Poissonova rozdělení: pravděpodobnost nejvýše k událostí se střední hodnotou λ je Q(k + 1, λ). Směrování výpočtu přes horní neúplnou funkci gamma namísto sčítání k + 1 Poissonových členů je opět volbou vyhodnocení jednoho konvergentního výrazu namísto kumulace mnoha malých hodnot
Diskrétní pravděpodobnosti bez přetečení faktoriálu
Diskrétní rozdělení přinášejí jiné riziko. Binomická pravděpodobnostní funkce (PMF) zahrnuje binomický koeficient, přičemž koeficient pro kombinaci dvacet šest z padesáti dvou je obrovské číslo. Pokud jej vytvoříte přímo, čitatel přeteče limit double dříve, než dělení vrátí výsledek zpět do smysluplného rozsahu pravděpodobnosti. Jádro jej proto přímo nevytváří. Počítá faktoriály v logaritmickém prostoru pomocí logaritmické funkce gamma, sčítá a odčítá logaritmy, zahrnuje logaritmus pravděpodobností úspěchu a neúspěchu a exponenciální funkci aplikuje až jednou na úplném konci
// Binomial probability mass, evaluated entirely in log space.
// LnGammaF(n+1) is ln(n!); the three log-factorials form ln(C(n,k)),
// and the whole exponent is built before a single Exp call.
// ln P(X=k) = ln(n!) - ln(k!) - ln((n-k)!) + k*ln(p) + (n-k)*ln(1-p)
result := Exp(LnGammaF(nt + 1) - LnGammaF(kk + 1) - LnGammaF(nt - kk + 1)
+ kk * Ln(pp) + (nt - kk) * Ln(1 - pp));Samotná logaritmická funkce gamma je Lanczosovou aproximací, přesnou na celé kladné ose a výpočetně nenáročnou. Protože každá velká veličina je držena jako svůj logaritmus až do závěrečného Exp, největším číslem, které procedura reálně vytvoří, je samotná pravděpodobnost, která je maximálně jedna. Poissonova hmotnostní funkce se řídí stejným receptem, přičemž jediný logaritmický člen gamma zastupuje faktoriál ve jmenovateli. Uzavřené tvary jsou speciálně ošetřeny na okrajích, kde je p přesně nula nebo jedna, takže kód nikdy nevolá Ln(0). HotXLS vrací hodnotu 0.2460938 pro BINOM.DIST(5,10,0.5,FALSE) a 0.6766764 pro kumulativní POISSON.DIST(2,2,TRUE), což odpovídá Excelu ve všech zobrazených číslicích
Inverzní funkce ohraničením dopředné křivky
Inverzní distribuční funkce se ptá na opačnou otázku: pro danou pravděpodobnost najděte hodnotu, jejíž CDF se jí rovná. Pouze jedna inverzní funkce v této sadě má rychlý přímý vzorec. NORM.S.INV, inverzní standardní normální rozdělení, používá Acklamovu racionální aproximaci, dvojici polynomických poměrů přesných přibližně na přesnost typu double v celém rozsahu, rozdělených na centrální oblast a dvě okrajové části. Jedná se o vyhodnocení v uzavřené formě bez iterace
Ostatní inverzní funkce takový vzorec nemají, a proto je jádro invertuje numericky. Ohraničí výsledek spodní a horní mezí zvolenou z definičního oboru rozdělení a poté provede bisekci (půlení intervalu): vyhodnotí dopřednou CDF ve středu intervalu, posune tu mez, která udržuje cílovou pravděpodobnost uzavřenou, a postup opakuje, dokud není interval dostatečně úzký. Pro inverzní funkce rozdělení gamma a chí-kvadrát začíná ohraničení na nule a velkorysém horním odhadu sestaveném z tvaru a měřítka, přičemž horní mez se zdvojnásobí, pokud pravděpodobnost ještě není uzavřena. Inverzní rozdělení t ohraničuje symetrické meze, které se rozšiřují směrem ven; inverzní rozdělení F provádí půlení na nezáporném intervalu. Náklady představují několik desítek vyhodnocení CDF na jedno volání, což je při rychlosti tabulkového procesoru neznatelné, a výhodou je, že každá inverzní funkce je přesně tak přesná, jako dopředná funkce, kterou invertuje. To je důvod, proč převod tam a zpět jako CHISQ.DIST(CHISQ.INV(0.7,5),5,TRUE) vrací hodnotu 0.7 s naprostou přesností
Dekadický logaritmus, který se skrýval v okrajové části
Acklamova procedura inverzního normálního rozdělení má tři větve. Široká centrální větev, používaná kdykoli se pravděpodobnost pohybuje přibližně mezi 0,025 a 0,975, zpracovává vstup pomocí polynomického poměru bez jakéhokoli logaritmu. Obě okrajové větve pro velmi malé nebo velmi velké pravděpodobnosti nejprve berou logaritmus vstupu, protože chování v okrajové části odpovídá odmocnině ze záporného přirozeného logaritmu p
Raná verze okrajové větvě používala dekadický logaritmus (o základu 10) tam, kam patřil logaritmus přirozený. Tyto dva se liší konstantním faktorem přibližně 2,30, takže výsledky v okrajových částech byly trvale chybně posunuté o značný rozdíl. Přesto funkce při běžné kontrole vypadala v pořádku, protože běžné kontroly se pohybují ve středu hodnot. NORM.S.INV(0.5) je nula, NORM.S.INV(0.975) je učebnicových 1.959964 a obě tyto hodnoty procházejí centrálním polynomem, který logaritmus vůbec nevolá. Chyba se projevila až v okamžiku, kdy pravděpodobnost přešla do okrajové části, například NORM.S.INV(0.001), což má vrátit -3.0902323 a namísto toho vracelo hodnotu zkreslenou o poměr přirozeného a dekadického logaritmu. Jakákoli funkce závislá na inverzním normálním rozdělení ve své okrajové části, včetně pomocných funkcí pro intervaly spolehlivosti, toto zkreslení zdědila. Poučení je jednoduché, ale cenné: funkce s větvenou strukturou vyžaduje testovací body uvnitř každé větve, protože správná běžná cesta snadno maskuje nefunkční ojedinělou větev. Nápravou byla změna jediného tokenu z dekadického logaritmu na přirozený logaritmus, čímž se okrajové hodnoty ihned srovnaly s hodnotami z Excelu
Znaménko x určuje okrajovou část t-rozdělení
Kumulativní funkce Studentova rozdělení t v sobě nese jemnost, kterou lze snadno obrátit. Její hodnota pochází z regularizované neúplné funkce beta vyhodnocené v bodě df / (df + x²), ale tato hodnota beta představuje pravděpodobnost v okrajové části mimo rozsah x a nikoli kumulativní pravděpodobnost do hodnoty x. Symetrický tvar t-rozdělení znamená, že převod závisí na tom, na jaké straně od nuly se x nachází
Pro x větší než nula je kumulativní pravděpodobnost jedna minus polovina symetrické okrajové části; pro x menší než nula je to polovina této okrajové části; v nule je to přesně polovina. Pokud byste hodnotu beta vrátili přímo, nahlásili byste špatnou stranu rozdělení, která se pro jakékoli nenulové x liší o celou plochu křivky. Pravostranná a dvoustranná varianta staví na stejné větvi, což je důvod, proč T.DIST.2T(1,1) vrací 0.5 a T.DIST(1,1,TRUE) vrací 0.75, a inverzní funkce T.INV pak provádí půlení intervalu vůči této opravené CDF, takže se převod tam i zpět uzavírá
Nic z toho není z buňky viditelné, což je také záměrem. Zapíšete vzorec a přečtete číslo, které odpovídá Excelu. Pokud rozšiřujete jádro o vlastní logiku, mechanika registrace funkce je popsána v našem průvodci jádrem vzorců a vlastními funkcemi, a způsob, jakým vzorce přistupují k ostatním listům a pojmenovaným rozsahům, popisuje článek o definovaných názvech a vzorcích napříč listy. Vše se dodává jako součást produktu HotXLS spreadsheet component pro Delphi a C++Builder, společně s rozhraními API pro čtení, zápis, grafy a formátování popsanými na jiných místech tohoto blogu