Írja be az =NORM.DIST(115,100,15,TRUE) képletet egy cellába, és az Excel különösebb felhajtás nélkül a 0.8413447 értéket adja vissza. A hívás egyszerű keresésnek tűnik. De nem az. Ezen egyetlen szám mögött a kumulatív normális eloszlás áll, egy zárt alakban nem integrálható függvény, a CHISQ.INV.RT és a BETA.DIST mögött pedig olyan speciális függvények rejtőznek, amelyeket egy gondos könyvtárnak ki kell számítania, nem pedig kézzel közelítenie. Egy Excel-kompatibilitást ígérő táblázatkezelő komponensnek ezeket az értékeket az Excel által mutatott utolsó számjegyig reprodukálnia kell, ami a numerikus módszerek, nem pedig csak a függvénynevek leképezését jelenti.
A HotXLS több mint ötvenet valósít meg ezek közül a statisztikai függvények közül, és a helyességüket biztosító munka szinte teljesen láthatatlan a képletsorból. Ez egy bemutató arról, hogyan számítja ki őket a motor: a megosztott speciális függvények magja, az aritmetika stabilitását megőrző elágazási döntések, valamint egy inverz-normális hiba, amely sokáig rejtve maradt a szélén (tail), mert a gyakori esetek soha nem érintették a hibás sort.
Egyetlen munkalaphívás, ötven eloszlás mögötte
A függvények átfogják azokat a családokat, amelyekért egy statisztikai munkafüzet nyúl. Ott van a normális eloszlás családja, a NORM.DIST és a NORM.S.DIST az inverzeikkel; a gamma és khí-négyzet család, a GAMMA.DIST, CHISQ.DIST, CHISQ.DIST.RT, CHISQ.INV.RT; a béta család, a BETA.DIST és a BETA.INV; a mintavételi eloszlások, mint a T.DIST, T.DIST.2T, F.DIST és F.INV; a diszkrét páros, a BINOM.DIST és a POISSON.DIST; valamint a következtetési segédek, mint a CONFIDENCE.T és a CONFIDENCE.NORM. A hívó szemszögéből mindegyik egyetlen képlet. Beállítja a bemeneteket a cellákban, kéri a munkafüzetet a kiértékelésre, és leolvassa az eredményt.
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;A munkafüzet Calculate metódusa lefordítja és kiértékeli az ad-hoc képletet az élő lapra nézve, és egy Variant-ot ad vissza. Egy részlet megzavarhatja a fejlesztőket az első próbálkozásnál: a Calculate mögötti képletértelmező a pontosvesszőt tekinti az argumentum elválasztójának, így a forma =SUM(A1;B1), nem pedig =SUM(A1,B1). A tárolt cellaképletek megtartják az Excel-szabványos vesszőt. Ugyanaz a kiértékelő indítja el az összes statisztikai függvényt az alábbiakban, így amint az egyik működik a Calculate-ben, a többi ugyanazt az utat követi.
A két függvény, amelyre minden más épül
A halmazban lévő kumulatív eloszlások többségét nem a saját definíciójuk összegzésével vagy integrálásával számítjuk ki. Két speciális függvényből számítjuk ki őket: a regularizált alsó inkomplett gammából (P(a, x)) és a regularizált inkomplett bétából (Ix(a, b)). Belsőleg ezek azok a segédfüggvények, amelyekre a diszpécserek támaszkodnak, és a lánc meglehetősen rövid. A khí-négyzet CDF a gamma CDF df/2 alakkal és 2-es skálával. A gamma CDF közvetlenül a P(a, x). A t, F és binomiális kumulatív függvények mind a regularizált inkomplett béta értékei a megfelelő argumentumoknál. A Poisson CDF a felső inkomplett gamma Q. Ha jól implementálja a gamma és béta függvényeket, egy tucat eloszlás örökli meg a pontosságukat ingyen.
A „regularizált” szó a lényeg. A nyers inkomplett gamma úgy növekszik, mint egy faktoriális, a nyers béta integrál pedig alul- vagy túlcsordulhat jóval azelőtt, hogy a válasz meglenne. A regularizált formák el vannak osztva a teljes gammával vagy bétával, így teljesen a nulla és egy közötti intervallumban élnek, ami pontosan az a tartomány, amelyet egy valószínűség elfoglal. Ez a normalizálás teszi lehetővé, hogy ugyanaz a rutin kiszolgáljon egy két szabadsági fokú khí-négyzetet és egy kétszáz szabadsági fokút anélkül, hogy a köztes tagok lefutnának a double végéről. Ez azt is megmagyarázza, miért nem számítunk ki CDF-et a sűrűség-tagok hosszú sorának összeadásával: minden tag hordozza a saját kerekítési hibáját, a hibák felhalmozódnak a sorozat futásakor, a regularizált speciális függvény pedig teljesen megkerüli az összeget egy gyorsan konvergáló sorozat vagy lánctört kiértékelésével.
Sorozat az átló alatt, lánctört felette
Az inkomplett gamma rutin egy döntést hoz, mielőtt bármit is kiszámolna: összehasonlítja az x-et az a + 1 értékkel. Ez a határ nem önkényes. A P(a, x) hatványsoros kifejtése gyorsan konvergál, ha az x kicsi az a-hoz képest, és lassan – végül használhatatlanul –, ha az x nagy. A lánctört ennek az ellenkezőjét mutatja. Így a motor a hatványsort használja az a + 1 alatti x értékekre, és Lentz lánctörtet az a + 1 feletti vagy egyenlő x értékekre, és mindegyik ágtól csak a saját feladatát kérjük el.
A lánctörtnek egy védelemre van szüksége. A Lentz-módszer úgy működik, hogy futó számlálót és nevezőt visz végig, és minden lépésben megfordítja a nevezőt, és ha bármelyik megközelíti a nullát, az inverzió felrobban. A javítás egy apró alsó korlát (floor): valahányszor egy köztes tag nagyságrendje nagyjából 1e-30 alá esik, rögzítjük 1e-30-ra, ami véges szinten tartja a rekurziót a konvergált érték megzavarása nélkül. Ugyanez a rögzítés jelenik meg az inkomplett béta lánctörtjében is ugyanezen okból. Ez egy kis állandó, amely komoly munkát végez, a különbséget jelentve a stabil kiértékelés és a nullától megkülönböztethetetlen számmal való osztás között.
A felső szél (upper tail), azaz a Q(a, x), egyszerűen 1 mínusz P(a, x), és a Poisson kumulatív ága így számítódik ki: a legfeljebb k esemény valószínűsége λ átlaggal a Q(k + 1, λ). Ennek átvezetése a felső inkomplett gammán ahelyett, hogy k + 1 Poisson tagot összegeznénk, ismét csak egy konvergens kifejezés kiértékelését jelenti sok kicsi felhalmozása helyett.
Diszkrét tömegek faktoriális túlcsordulás nélkül
A diszkrét eloszlások más veszélyt jelentenek. A binomiális valószínűségi tömeg (probability mass) binomiális együtthatót tartalmaz, és az 52-ből 26 kiválasztásának együtthatója hatalmas egész szám. Ha közvetlenül hozzuk létre, a számláló túlcsordul a double-ben az osztás előtt, amely visszahozná azt egy ésszerű valószínűségre. A motor soha nem hozza létre közvetlenül. A faktoriálisokat logaritmikus térben számítja ki a log-gamma függvény segítségével, összeadja és kivonja a logaritmusokat, bevonja a siker- és kudarcvalószínűségek logaritmusát, és a legvégén egyszer hatványoz (exponentiates).
// 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));Maga a log-gamma függvény egy Lánczos-féle közelítés, amely a teljes pozitív tengelyen pontos és olcsón értékelhető ki. Mivel minden nagy mennyiség a logaritmusaként van tárolva a végső Exp hívásig, a legnagyobb szám, amelyet a rutin valaha előállít, maga a valószínűség, amely legfeljebb egy. A Poisson tömegfüggvény ugyanezt a receptet követi, az egyetlen log-gamma taggal helyettesítve a faktoriálist a nevezőben. A zárt formák speciális esetként kezeltek a határokon, ahol a p pontosan nulla vagy egy, így a kód soha nem hívja a Ln(0)-t. A HotXLS 0.2460938-at ad vissza a BINOM.DIST(5,10,0.5,FALSE) függvényre és 0.6766764-et a kumulatív POISSON.DIST(2,2,TRUE)-ra, egyezve az Excel által kiírt számjegyekkel.
Inverzek az előre görbe behatárolásával
Az inverz eloszlásfüggvény az ellenkező kérdést teszi fel: adott valószínűséghez keresse meg azt az értéket, amelynek a CDF-je egyenlő vele. Ebben a halmazban csak egyetlen inverz rendelkezik gyors, közvetlen képlettel. A NORM.S.INV, az inverz standard normális eloszlás, egy Acklam racionális közelítést használ, amely egy polinomarány-pár, amely nagyjából egy double pontosságával rendelkezik a teljes tartományban, egy középső régióra és két szélre bontva. Ez egy zárt alakú kiértékelés iteráció nélkül.
A többi inverz nem rendelkezik ilyen képlettel, így a motor numerikusan invertálja őket. Közrefogja (bracket) a választ az eloszlás tartományából választott alsó és felső határral, majd felez (bisects): kiértékeli az előre CDF-et a középpontban, elmozdítja azt a határt, amely a megcélzott valószínűséget közrezárva tartja, és addig ismétli, amíg az intervallum szűk nem lesz. A gamma és khí-négyzet inverzeknél a közrefogás nullánál kezdődik, és az alakból és skálából épített bőséges felső becslésnél, megduplázva a felső határt, ha a valószínűség még nincs közrefogva. A t-inverz szimmetrikus határokat fog közre, amelyek kifelé tágulnak; az F-inverz nemnegatív intervallumon felez. A költség hívásonként néhány tucat CDF-kiértékelés, ami észrevehetetlen a táblázatkezelés sebességével, az előny pedig az, hogy minden inverz pontosan olyan pontos, mint az előre függvény, amelyet invertál. Ezért ad vissza egy olyan körút, mint a CHISQ.DIST(CHISQ.INV(0.7,5),5,TRUE) hajszálpontosan 0.7-et.
A tízes alapú logaritmus, amely a szélben bujkált
Íme a hiba, amelyet érdemes elmesélni, mert az a fajta, amely sokáig fennmarad. Az Acklam inverz-normál rutinnak három ága van. A széles középső ág, amely akkor használatos, ha a valószínűség körülbelül 0.025 és 0.975 között van, a bemenetet egy polinomarányon futtatja át, amelyben nincs sehol logaritmus. A két szélső (tail) ág a nagyon kicsi vagy nagyon nagy valószínűségekhez először a bemenet logaritmusát veszi, mert a szélek úgy viselkednek, mint a p természetes logaritmusának mínusz egyszeresének négyzetgyöke.
A szélső ág egy korai verziója tízes alapú logaritmust vett ott, ahová a természetes logaritmus tartozott. A kettő egy körülbelül 2.30-as állandó szorzóval tér el egymástól, így a szélső eredmények konzisztens, jelentős mértékben hibásak voltak. A függvény mégis jónak tűnt minden alkalmi ellenőrzésnél, mert az alkalmi ellenőrzések középen zajlanak. A NORM.S.INV(0.5) értéke nulla, a NORM.S.INV(0.975) a tankönyvi 1.959964, és mindkettő azon a középső polinomon fut át, amely egyáltalán nem hív logaritmust. A hiba csak akkor jelentkezett, ha a valószínűség átlépett a szélre – például a NORM.S.INV(0.001) esetében, amelynek -3.0902323-at kell visszaadnia, ehelyett a természetes logaritmus kontra tízes alapú logaritmus eltérésével tért vissza. Bármely függvény, amely a széleken az inverz normálistól függ – beleértve a konfidenciaintervallum-segédeket is –, örökölte ugyanezt a torzulást. A tanulság hétköznapi és drága: egy elágazó struktúrájú függvénynek minden ágban szüksége van tesztpontokra, mert a helyes gyakori útvonal vidáman elfed egy hibás ritkát. A javítás egyetlen token változtatása volt a tízes alapú logaritmusról a természetes logaritmusra, és a szélső értékek azonnal az Excel értékeihez igazodtak.
Az x előjele dönti el a t-eloszlás szélét
A Student t kumulatív függvény olyan finomságot hordoz, amelyet könnyű visszafelé érteni. Az értéke a regularizált inkomplett bétából származik, amelyet a df / (df + x²) helyen értékelünk ki, de ez a béta érték az x nagyságán túli szél valószínűsége, nem pedig a kumulatív valószínűség x-ig. A t-eloszlás szimmetrikus alakja azt jelenti, hogy a konverzió attól függ, hogy az x a nulla melyik oldalára esik.
// 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 meanNulla feletti x esetén a kumulatív valószínűség egy mínusz a szimmetrikus szél fele; nulla alatti x esetén a szél fele; nullánál pedig pontosan egyketted. Ha közvetlenül a béta értéket adja vissza, az eloszlás rossz oldalát jelenti, a görbe teljes törzsével eltolódva bármely nemnulla x esetén. A jobb szélső és két szélső változatok ugyanerre az ágra épülnek, ezért tér vissza a T.DIST.2T(1,1) 0.5-ként, a T.DIST(1,1,TRUE) pedig 0.75-ként, és az inverz T.INV ezen korrigált CDF-fel szemben felez, így zárul a körút.
Mindebből semmi sem látható a cellából, és ez a tervezett eredmény. Ön beír egy képletet, és leolvas egy számot, amely megegyezik az Excellel. Ha saját logikával egészíti ki a motort, a függvény regisztrálásának mechanikáját a képletmotorról és egyedi függvényekről szóló bemutatónk ismerteti, a képletek lapok és elnevezett tartományok közötti elérésének módját pedig a definiált nevekről és keresztlap-képletekről szóló cikk tárgyalja. Mindez a Delphi és C++Builder platformokhoz készült HotXLS spreadsheet component részeként érhető el, a blogunkon máshol ismertetett olvasási, írási, diagramkészítési és formázási API-kkal együtt.