Zadejte do buňky =NORM.DIST(115,100,15,TRUE) a Excel bez okolků vrátí 0.8413447. Volání vypadá jako jednoduché vyhledání. 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 CHISQ.INV.RT a BETA.DIST stojí speciální funkce, které pečlivá knihovna musí skutečně vyhodnotit, nikoli jen odhadnout od ruky. Tabulkový komponent, který si nárokuje kompatibilitu s Excelem, musí tyto hodnoty reprodukovat do poslední číslice, kterou Excel zobrazuje, což znamená reprodukovat numerické metody, nikoli jen názvy funkcí.
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ěř nepostřehnutelná. Toto je přehled toho, jak je engine počítá: sdílené jádro speciálních funkcí, rozhodování o větvení udržující stabilní aritmetiku a jedna chyba inverzního normálního rozdělení, která se dlouho skrývala v okrajových částech (tail), protože běžné případy se na poškozený řádek kódu nikdy nedostaly.
Jedno volání v listu, padesát rozdělení za ním
Funkce pokrývají rodiny, po kterých statistický sešit sahá. Je zde normální rozdělení, NORM.DIST a NORM.S.DIST s jejich inverzemi; rozdělení gama 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íci pro statistickou indukci jako CONFIDENCE.T a CONFIDENCE.NORM. Z pohledu volajícího je každá funkce 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 živému listu a vrací Variant. Jeden detail lidi při prvním pokusu zaskočí: parser vzorců na pozadí metody Calculate používá jako oddělovač argumentů středník, takže se píše =SUM(A1;B1), nikoli =SUM(A1,B1). Uložené vzorce buněk si ponechávají standardní čárku z Excelu. Stejný vyhodnocovač odbavuje každou níže uvedenou statistickou funkci, 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 distribučních funkcí v této sadě se nepočítá sčítáním nebo integrací jejich vlastních definic. Počítají se ze dvou speciálních funkcí: regularizované dolní neúplné gama funkce, zapisované jako P(a, x), a regularizované neúplné beta funkce, zapisované jako Ix(a, b). Interně se jedná o pomocníky, o které se vyhodnocovací moduly opírají, a cesta k nim je krátká. CDF (kumulativní distribuční funkce) chí-kvadrát je CDF gama s tvarem df/2 a měřítkem 2. CDF gama je přímo P(a, x). Kumulativní funkce rozdělení t, F a binomického jsou všechny hodnotami regularizované neúplné beta funkce při správných argumentech. CDF Poissonova rozdělení je horní neúplná gama Q. Implementujte dobře funkce gama a beta a tucet rozdělení zdědí jejich přesnost zdarma.
Slovo "regularizovaná" je tím nejdůležitějším. Původní neúplná gama roste jako faktoriál a původní integrál beta může podtéct nebo přetéct dlouho předtím, než se dopracujete k výsledku. Regularizované formy jsou vyděleny kompletní funkcí gama nebo beta, takže žijí zcela v intervalu od nuly do jedné, což je přesně rozsah, který zaujímá pravděpodobnost. Tato normalizace umožňuje, aby stejná rutina sloužila pro chí-kvadrát se dvěma stupni volnosti i se dvěma sty, aniž by mezilehlé členy přetekly mimo rozsah double. Vysvětluje to také, proč nepočítáte CDF sčítáním dlouhé řady hustotních členů: každý člen nese svou vlastní zaokrouhlovací chybu, chyby se při běhu řady sčítají a regularizovaná speciální funkce se součtu zcela vyhne vyhodnocením rychle konvergující řady nebo řetězového zlomku.
Řada pod diagonálou, řetězový zlomek nad ní
Rutina neúplné gama funkce dělá před jakýmkoli výpočtem jedno rozhodnutí: porovnává x s hodnotou a + 1. Tato hranice není libovolná. Rozvoj mocninné řady P(a, x) konverguje rychle, když je x malé vzhledem k a, a pomalu (nakonec nepoužitelně), když je x velké. Řetězový zlomek má opačný charakter. Engine 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, takže každá větev provádí pouze práci, pro kterou je vhodná.
Řetězový zlomek vyžaduje jednu ochranu. Lentzova metoda funguje tak, že nese průběžný čitatel a jmenovatel a v každém kroku jmenovatel invertuje; pokud se kterýkoli z nich přiblíží nule, inverze selže. Nápravou je drobná dolní mez (floor): kdykoli velikost mezilehlého členu klesne pod přibližně 1e-30, je omezena na hodnotu 1e-30, což udržuje rekurenci konečnou, aniž by to narušilo konvergovanou hodnotu. Stejné omezení se ze stejného důvodu objevuje v řetězovém zlomku neúplné beta funkce. Jedná se o malou konstantu plnící důležitou roli, což představuje rozdíl mezi stabilním vyhodnocením a dělením hodnotou nerozeznatelnou od nuly.
Horní konec (upper tail), Q(a, x), je jednoduše 1 minus P(a, x), a tak se počítá kumulativní větev Poissonova rozdělení: pravděpodobnost nejvýše k událostí se střední hodnotou λ je Q(k + 1, λ). Jeho nasměrování přes horní neúplnou gama funkci namísto sčítání k + 1 Poissonových členů je opět volbou vyhodnotit jeden konvergentní výraz namísto akumulování mnoha malých hodnot.
Diskrétní rozdělení bez přetečení faktoriálu
Diskrétní rozdělení přinášejí jiné riziko. Binomická pravděpodobnostní funkce (probability mass) zahrnuje binomický koeficient a koeficient pro kombinaci 52 nad 26 je obrovské celé číslo. Pokud jej vytvoříte přímo, čitatel přeteče rozsah double ještě před dělením, které by jej vrátilo zpět na rozumnou pravděpodobnost. Engine jej nikdy nevytváří přímo. Počítá faktoriály v logaritmickém prostoru pomocí logaritmické gama funkce, sčítá a odečítá logaritmy, přičítá logaritmus pravděpodobností úspěchu a neúspěchu a exponenciální funkci aplikuje až na samé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á gama funkce je Lanczosovou aproximací, přesnou na celé kladné ose a levnou na výpočet. 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é rutina kdy zhmotní, je samotná pravděpodobnost, která je nejvýše jedna. Poissonova pravděpodobnostní funkce sleduje stejný recept, přičemž jediný člen logaritmické gama zastupuje faktoriál ve jmenovateli. Uzavřené formy jsou ošetřeny jako speciální případy na okrajích, kde je p přesně nula nebo jedna, takže kód nikdy nevolá Ln(0). HotXLS vrací 0.2460938 pro BINOM.DIST(5,10,0.5,FALSE) a 0.6766764 pro kumulativní POISSON.DIST(2,2,TRUE), což odpovídá hodnotám v Excelu v celém rozsahu zobrazených číslic.
Inverze ohraničením přímé křivky
Inverzní distribuční funkce klade opačnou otázku: pro danou pravděpodobnost najděte hodnotu, jejíž CDF se jí rovná. Pouze jedna inverze 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 podílů přesných zhruba na přesnost double v celém rozsahu), rozdělenou na centrální oblast a dva okraje (tails). Jde o vyhodnocení v uzavřeném tvaru bez iterace.
Ostatní inverze takový vzorec nemají, takže je engine invertuje numericky. Ohraničí (brackets) odpověď dolní a horní mezí zvolenou z definičního oboru rozdělení a poté provádí bisekci (půlení intervalu): vyhodnotí přímou CDF ve středu, posune tu mez, která udržuje cílovou pravděpodobnost uzavřenou, a opakuje to, dokud není interval úzký. U inverzí gama 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. Inverze t ohraničuje symetrické meze, které se rozšiřují směrem ven; inverze F provádí bisekci na nezáporném intervalu. Nákladem je několik desítek vyhodnocení CDF na volání, což je při rychlosti tabulkového procesoru nepostřehnutelné, a výhodou je, že každá inverze je přesně tak přesná, jako přímá funkce, kterou invertuje. Proto obousměrný převod jako CHISQ.DIST(CHISQ.INV(0.7,5),5,TRUE) vrací 0.7 s vysokou přesností.
Dekadický logaritmus, který se skrýval v okraji
Zde je chyba, kterou stojí za to popsat, protože je to ten typ, který přežívá dlouho. Acklamova rutina pro inverzní normální rozdělení má tři větve. Široká centrální větev, používaná, kdykoli pravděpodobnost leží přibližně mezi 0.025 a 0.975, zpracovává vstup přes polynomický podíl bez jakéhokoli logaritmu. Obě okrajové větve, pro velmi malé nebo velmi velké pravděpodobnosti, nejprve přebírají logaritmus vstupu, protože okraj se chová jako druhá odmocnina z minus přirozeného logaritmu p.
Raná verze okrajové větve používala dekadický logaritmus tam, kam patřil přirozený logaritmus. Tyto dva se liší konstantním faktorem přibližně 2.30, so okrajové výsledky byly wrong o konzistentní, znatelný kus. Přesto funkce při každé běžné kontrole vypadala v pořádku, protože běžné kontroly se provádějí uprostřed. 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 okraje, například NORM.S.INV(0.001), což musí vrátit -3.0902323 a místo toho se vrátila hodnota posunutá 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ém okraji, včetně pomocníků pro intervaly spolehlivosti, zdědila stejné zkreslení. Ponaučení je všední a drahé: funkce s větvenou strukturou potřebuje testovací body v každé větvi, protože správná běžná cesta snadno maskuje nefunkční vzácnou větev. Nápravou byla změna jednoho tokenu z dekadického logaritmu na přirozený a hodnoty na okrajích se okamžitě srovnaly s Excelem.
Znaménko x rozhoduje o okraji t-rozdělení
Kumulativní funkce Studentova t-rozdělení nese záludnost, kterou lze snadno obrátit. Její hodnota pochází z regularizované neúplné beta funkce vyhodnocené při df / (df + x²), ale tato hodnota beta je pravděpodobností v okraji za hodnotou x, nikoli kumulativní pravděpodobností do x. Symetrický tvar t-rozdělení znamená, že převod závisí na tom, na které straně od nuly se x nachází.
// Student t CDF. ib is the regularized incomplete beta at df/(df+x*x),
// which measures the symmetric tail. The cumulative value depends on
// the sign of x; returning ib unconverted gives the wrong tail.
ib := BetaIF(df / 2, 0.5, df / (df + x * x));
if x > 0 then
result := 1 - 0.5 * ib // above the mean: one minus half the tail
else if x < 0 then
result := 0.5 * ib // below the mean: half the tail
else
result := 0.5; // exactly at the meanPro x větší než nula je kumulativní pravděpodobnost rovna jedna minus polovina symetrického okraje; pro x menší než nula je to polovina tohoto okraje; na nule je to přesně jedna polovina. Pokud vrátíte hodnotu beta přímo, ohlásíte špatnou stranu rozdělení, což je pro jakékoli nenulové x posun o celé tělo křivky. Varianty pro pravý okraj a dva okraje staví na stejné větvi, což je důvod, proč se T.DIST.2T(1,1) vrací jako 0.5 a T.DIST(1,1,TRUE) jako 0.75 a inverzní T.INV provádí bisekci vůči této opravené CDF, čímž se obousměrný převod uzavírá.
Nic z toho není z buňky viditelné, což je zamýšlený výsledek. Zapíšete vzorec a přečtete číslo, které souhlasí s Excelem. Pokud rozšiřujete engine o vlastní logiku, mechanika registrace funkce je popsána v našem průvodci formulovým enginem a vlastními funkcemi a způsob, jakým vzorce přistupují k listům a pojmenovaným rozsahům, rozebírá článek o definovaných názvech a vzorcích napříč listy. Vše je dodáváno v produktu HotXLS tabulkový komponent pro Delphi a C++Builder spolu s rozhraními API pro čtení, zápis, grafy a formátování popsanými na jiných místech tohoto blogu.