Zapíšte =NORM.DIST(115,100,15,TRUE) do bunky a Excel bez rečí vráti 0.8413447. Volanie vyzerá ako jednoduché vyhľadanie. Nie je to tak. Za týmto jediným číslom je kumulatívne normálne rozdelenie, integrál bez uzavretého tvaru, a za CHISQ.INV.RT a BETA.DIST sedia špeciálne funkcie, ktoré musí starostlivá knižnica vyhodnotiť, a nie približne odhadovať ručne. Tabuľkový komponent, ktorý si nárokuje kompatibilitu s Excelom, musí tieto hodnoty reprodukovať do poslednej číslice, ktorú Excel zobrazuje, čo znamená reprodukovanie numerických metód, nielen názvov funkcií.
HotXLS implementuje viac ako päťdesiat týchto štatistických funkcií a práca, ktorá ich robí správnymi, je z riadku vzorcov takmer úplne neviditeľná. Toto je prehliadka toho, ako ich engine počíta: zdieľané jadro špeciálnych funkcií, rozhodnutia o vetvení, ktoré udržujú aritmetiku stabilnú, a jedna chyba inverzného normálneho rozdelenia, ktorá sa dlho skrývala na okraji (chvoste), pretože bežný prípad sa nefunkčného riadku nikdy nedotkol.
Jedno volanie v hárku, päťdesiat rozdelení za ním
Funkcie pokrývajú rodiny, po ktorých štatistický zošit siaha. Je tu rodina normálneho rozdelenia, NORM.DIST a NORM.S.DIST s ich inverziami; rodina gama a chí-kvadrát rozdelenia, GAMMA.DIST, CHISQ.DIST, CHISQ.DIST.RT, CHISQ.INV.RT; rodina beta rozdelenia, BETA.DIST a BETA.INV; výberové rozdelenia T.DIST, T.DIST.2T, F.DIST a F.INV; diskrétna dvojica BINOM.DIST and POISSON.DIST; a pomocníci pre odvodzovanie ako CONFIDENCE.T a CONFIDENCE.NORM. Z pohľadu volajúceho je každá z nich jediným vzorcom. Nastavíte vstupy do buniek, požiadate zošit o vyhodnotenie a prečítate výsledok.
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;Metóda Calculate zošita skompiluje a vyhodnotí ad-hoc vzorec voči živému hárku a vráti Variant. Jeden detail ľudí pri prvom pokuse prekvapí: parser vzorcov za metódou Calculate používa ako oddeľovač argumentov bodkočiarku, takže ide o =SUM(A1;B1), nie =SUM(A1,B1). Uložené vzorce buniek si zachovávajú štandardnú excelovskú čiarku. Rovnaký vyhodnocovač odosiela každú štatistickú funkciu nižšie, takže akonáhle jedna z nich funguje v Calculate, ostatné nasledujú rovnakú cestu.
Dve funkcie, na ktorých stojí všetko ostatné
Väčšina kumulatívnych distribúcií v tejto sade sa nepočíta sčítavaním alebo integrovaním ich vlastných definícií. Počítajú sa z dvoch špeciálnych funkcií: regularizovanej dolnej neúplnej gama funkcie, zapísanej ako P(a, x), a regularizovanej neúplnej beta funkcie, zapísanej ako Ix(a, b). Interne ide o pomocníkov, o ktorých sa dispečeri opierajú, a reťazec je krátky. CDF chí-kvadrát rozdelenia je CDF gama rozdelenia s tvarom df/2 a mierkou 2. CDF gama rozdelenia je priamo P(a, x). Kumulatívne funkcie rozdelení t, F a binomického sú všetky hodnotami regularizovanej neúplnej bety pri správnych argumentoch. Poissonova CDF je horná neúplná gama Q. Ak implementujete funkcie gama a beta správne, tucet rozdelení zdedí ich presnosť zadarmo.
Slovo "regularizovaný" je celým zmyslom. Surová neúplná gama rastie ako faktoriál a surový integrál bety môže podtiecť alebo pretiecť oveľa skôr, ako výsledok. Regularizované formy sú vydelené kompletnou gamou alebo betou, takže žijú výlučne v intervale od nuly po jeden, čo je presne rozsah, ktorý zaberá pravdepodobnosť. Táto normalizácia umožňuje, aby rovnaká rutina slúžila chí-kvadrátu s dvoma stupňami voľnosti aj s dvoma stovkami bez toho, aby stredné členy utiekli z rozsahu typu double. Vysvetľuje to tiež, prečo nepočítate CDF sčítaním dlhého chvosta hustotných členov: každý člen nesie svoju vlastnú zaokrúhľovaciu chybu, chyby sa hromadia pri behu radu a regularizovaná špeciálna funkcia súčet úplne obchádza vyhodnotením rýchlo konvergujúceho radu alebo reťazového zlomku.
Rad pod diagonálou, reťazový zlomok nad ňou
Rutina neúplnej gamy urobí jedno rozhodnutie pred akýmkoľvek výpočtom: porovná x s hodnotou a + 1. Táto hranica nie je ľubovoľná. Rozvoj P(a, x) do mocninového radu konverguje rýchlo, keď je x malé vzhľadom na a, a pomaly, prípadne nepoužiteľne, keď je x veľké. Reťazový zlomok má opačný charakter. Engine preto používa mocninový rad pre x menšie ako a + 1 a Lentzov reťazový zlomok pre x rovné alebo väčšie ako a + 1, a každá vetva je požiadaná iba o prácu, v ktorej je efektívna.
Reťazový zlomok potrebuje jednu ochranu. Lentzova metóda funguje tak, že prenáša priebežný čitateľ a menovateľ a v každom kroku invertuje menovateľ, a ak sa niektorý z nich priblíži k nule, inverzia zlyhá. Riešením je drobný spodný limit (floor): kedykoľvek stredný člen klesne pod veľkosť približne 1e-30, upne sa na 1e-30, čo udržuje rekurziu konečnú bez narušenia konvergovanej hodnoty. Rovnaké upnutie sa z rovnakého dôvodu objavuje v reťazovom zlomku neúplnej bety. Je to malá konštanta vykonávajúca dôležitú prácu, rozdiel medzi stabilným vyhodnotením a delením niečím, čo je nerozoznateľné od nuly.
Horný chvost, Q(a, x), je jednoducho 1 mínus P(a, x), a takto sa počíta Poissonova kumulatívna vetva: pravdepodobnosť najviac k udalostí so strednou hodnotou λ je Q(k + 1, λ). Jej smerovanie cez hornú neúplnú gamu namiesto sčítavania k + 1 Poissonových členov je opäť voľbou vyhodnotiť jeden konvergentný výraz namiesto akumulovania mnohých malých.
Diskrétne hmotnosti bez pretečenia faktoriálu
Diskrétne rozdelenia prinášajú iné riziko. Binomická pravdepodobnosť zahŕňa binomický koeficient a koeficient pre kombinačné číslo päťdesiatdva nad dvadsaťšesť je obrovské celé číslo. Ak ho vytvoríte priamo, čitateľ pretečie cez rozsah double ešte pred delením, ktoré by ho prinieslo späť na rozumnú pravdepodobnosť. Engine ho nikdy nevytvára priamo. Počíta faktoriály v logaritmickom priestore pomocou funkcie log-gama, sčítava a odčítava logaritmy, pripočítava logaritmus pravdepodobností úspechu a zlyhania a umocňuje až na samom 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á funkcia log-gama je Lanczosova aproximácia, presná na celej kladnej osi a lacná na vyhodnotenie. Keďže každá veľká veličina sa uchováva ako jej logaritmus až po finálne Exp, najväčšie číslo, ktoré rutina kedy materializuje, je samotná pravdepodobnosť, ktorá je najviac jedna. Poissonova hmotnostná funkcia sa riadi rovnakým receptom, pričom jeden člen log-gama zastupuje faktoriál v menovateli. Uzavreté tvary sú špeciálne ošetrené na okrajoch, kde p je presne nula alebo jedna, takže kód nikdy nevolá Ln(0). HotXLS vracia 0.2460938 pre BINOM.DIST(5,10,0.5,FALSE) a 0.6766764 pre kumulatívne POISSON.DIST(2,2,TRUE), čo sa zhoduje s Excelom vo všetkých zobrazených čísliciach.
Inverzné hodnoty ohraničením priamej krivky
Inverzná funkcia rozdelenia kladie opačnú otázku: pre danú pravdepodobnosť nájdite hodnotu, ktorej CDF sa jej rovná. Iba jedna inverzia v tejto sade má rýchly priamy vzorec. NORM.S.INV, inverzné štandardné normálne rozdelenie, používa Acklamovu racionálnu aproximáciu, dvojicu polynomiálnych pomerov presných približne na presnosť typu double v celom rozsahu, rozdelenú na strednú oblasť a dva chvosty. Ide o vyhodnotenie v uzavretom tvare bez iterácie.
Ostatné inverzie takýto vzorec nemajú, preto ich engine invertuje numericky. Ohraničí odpoveď spodnou a hornou hranicou vybranou z nosiča rozdelenia a potom robí bisekciu (polenie intervalu): vyhodnotí priamu CDF v strede, posunie tú hranicu, ktorá udržuje cieľovú pravdepodobnosť uzavretú, a opakuje, kým nie je interval úzky. Pre inverzie gamy a chí-kvadrátu začína ohraničenie na nule a veľkorysom hornom odhade postavenom z tvaru a mierky, pričom hornú hranicu zdvojnásobí, ak pravdepodobnosť ešte nie je uzavretá. Inverzia t ohraničuje symetrické hranice, ktoré sa rozširujú smerom von; inverzia F polí interval na nezápornom intervale. Náklady predstavujú niekoľko desiatok vyhodnotení CDF na volanie, čo je pri rýchlosti tabuľkového kalkulátora nepostrehnuteľné, a výhodou je, že každá inverzia je presne taká presná ako priama funkcia, ktorú invertuje. Preto spiatočná cesta ako CHISQ.DIST(CHISQ.INV(0.7,5),5,TRUE) vráti 0.7 s presnosťou na vlas.
Desiatkový logaritmus, ktorý sa skrýval v chvoste
Tu je chyba, ktorú stojí za to spomenúť, pretože patrí k tým, ktoré prežívajú dlho. Acklamova rutina inverzného normálneho rozdelenia má tri vetvy. Široká stredná vetva, používaná vždy, keď pravdepodobnosť leží približne medzi 0.025 a 0.975, preháňa vstup cez polynomiálny pomer bez akéhokoľvek logaritmu. Dve chvostové vetvy pre veľmi malé alebo veľmi veľké pravdepodobnosti berú najprv logaritmus vstupu, pretože chvost sa správa ako druhá odmocnina zo záporného prirodzeného logaritmu p.
Skorá verzia chvostovej vetvy brala desiatkový logaritmus tam, kam patril prirodzený logaritmus. Tieto dva sa líšia konštantným faktorom približne 2.30, sože výsledky na chvostoch boli nesprávne o konzistentnú, značnú odchýlku. A predsa funkcia vyzerala v poriadku pri každej bežnej kontrole, pretože bežné kontroly žijú v strede. NORM.S.INV(0.5) je nula, NORM.S.INV(0.975) je učebnicových 1.959964 a obe tieto hodnoty prechádzajú cez stredový polynóm, ktorý logaritmus vôbec nevolá. Chyba sa objavila až vtedy, keď pravdepodobnosť prešla do chvosta, napríklad NORM.S.INV(0.001), ktoré musí vrátiť -3.0902323 a namiesto toho sa vrátilo s odchýlkou kvôli pomeru prirodzeného a desiatkového logaritmu. Akákoľvek funkcia, ktorá vo svojom chvoste závisí od inverzného normálneho rozdelenia, vrátane pomocníkov pre intervaly spoľahlivosti, zdedila rovnaké skreslenie. Ponaučenie je obyčajné a drahé: funkcia so štruktúrou vetvenia potrebuje testovacie body v každej vetve, pretože správna bežná cesta veselo zamaskuje nefunkčnú zriedkavú. Nápravou bola zmena jedného tokenu z desiatkového logaritmu na prirodzený a hodnoty na chvostoch okamžite skočili na excelovské.
Znamienko x rozhoduje o chvoste t-rozdelenia
Kumulatívna funkcia Studentovho t-rozdelenia nesie jemnosť, ktorú je ľahké pochopiť naopak. Jej hodnota pochádza z regularizovanej neúplnej bety vyhodnotenej v bode df / (df + x²), ale táto hodnota bety je pravdepodobnosťou v chvoste za veľkosťou x, nie kumulatívnou pravdepodobnosťou do x. Symetrický tvar t-rozdelenia znamená, že prevod závisí od toho, na ktorej strane nuly x leží.
// 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 meanPre x nad nulou je kumulatívna pravdepodobnosť jedna mínus polovica symetrického chvosta; pre x pod nulou je to polovica tohto chvosta; v nule je to presne jedna polovica. Vráťte hodnotu bety priamo a nahlásite nesprávnu stranu rozdelenia, posunutú o celé telo krivky pre akékoľvek nenulové x. Pravostranný variant a variant s dvoma chvostami stavajú na rovnakej vetve, preto sa T.DIST.2T(1,1) vracia ako 0.5 a T.DIST(1,1,TRUE) ako 0.75, pričom inverzná T.INV robí bisekciu voči tejto opravenej CDF, takže spiatočná cesta sa uzavrie.
Nič z toho nie je z bunky viditeľné a to je zamýšľaný výsledok. Napíšete vzorec a prečítate číslo, ktoré súhlasí s Excelom. Ak rozširujete engine o vlastnú logiku, mechanizmus registrácie funkcie je popísaný v našom sprievodcovi vzorcovým enginom a vlastnými funkciami a spôsob, akým vzorce siahajú naprieč hárkami a pomenovanými rozsahmi, je opísaný v článku o definovaných názvoch a vzorcoch naprieč hárkami. Všetky veci sa dodávajú v rámci tabuľkového komponentu HotXLS spreadsheet component pre Delphi a C++Builder spolu s rozhraniami API pre čítanie, zápis, grafy a formátovanie, o ktorých sa píše inde na tomto blogu.