Vpišite =NORM.DIST(115,100,15,TRUE) v celico in Excel bo brez oklevanja vrnil 0.8413447. Klic deluje kot preprosto iskanje. Pa ni. Za to eno številko se skriva kumulativna normalna porazdelitev, integral brez zaključene oblike, za CHISQ.INV.RT in BETA.DIST pa stojijo posebne funkcije, ki jih mora natančna knjižnica ovrednotiti in ne le ročno približati. Komponenta preglednic, ki trdi, da je združljiva z Excelom, mora te vrednosti reproducirati do zadnjega mesta, ki ga prikazuje Excel, kar pomeni reprodukcijo numeričnih metod in ne le imen funkcij.
HotXLS implementira več kot petdeset teh statističnih funkcij, delo, ki zagotavlja njihovo pravilnost, pa je skoraj v celoti nevidno iz vrstice za formule. Temu je namenjen pregled, kako jih pogon računa: skupno jedro posebnih funkcij, odločitve o vejah, ki ohranjajo aritmetiko stabilno, in ena napaka inverzne normale, ki se je dolgo skrivala v repu, ker se pogostejši primeri niso nikoli dotaknili napačne vrstice.
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 delovnem zvezku prevede in ovrednoti ad-hoc formulo glede na aktivni list in vrne Variant. Ena podrobnost lahko uporabnike zmede ob prvem poskusu: razčlenjevalnik formul za metodo Calculate uporablja podpičje kot ločilo argumentov, torej =SUM(A1;B1) in ne =SUM(A1,B1). Shranjene formule v celicah ohranjajo Excelov standard z vejico. Isti ocenjevalnik odpošlje vsako spodnjo statistično funkcijo, tako da, ko ena od njih deluje v metodi Calculate, ostale sledijo isti poti.
Dve funkciji, na katerih gradi vse ostalo
Večina kumulativnih porazdelitev v tem naboru se ne računa s seštevanjem ali integracijo lastnih definicij. Računajo se iz dveh posebnih funkcij: regularizirane spodnje nepopolne gama (regularized lower incomplete gamma), zapisane kot P(a, x), in regularizirane nepopolne beta (regularized incomplete beta), zapisane kot Ix(a, b). Interno sta to pomočnika, na katera se opirajo dispečerji, veriga pa je kratka. Kumulativna funkcija porazdelitve (CDF) hi-kvadrat je gama CDF z obliko df/2 in merilom 2. Gama CDF je neposredno P(a, x). Kumulativne funkcije t, F in binomska porazdelitev so vse vrednosti regularizirane nepopolne beta funkcije pri ustreznih argumentih. Poissonova CDF je zgornja nepopolna gama Q. Dobro implementirajte funkciji gama in beta in ducat porazdelitev bo brezplačno podedovalo njihovo natančnost.
Beseda "regularizirana" predstavlja bistvo. Surova nepopolna gama raste kot fakulteta, surov beta integral pa lahko doživi podprekoračitev ali prekoračitev precej pred rezultatom. Regularizirane oblike so deljene s celotno gama ali beta funkcijo, zato živijo v celoti v intervalu od nič do ena, kar je natanko obseg, ki ga zaseda verjetnost. Ta normalizacija omogoča, da ista rutina služi hi-kvadratu z dvema prostostnima stopnjama in tistemu z dvesto, ne da bi vmesni členi presegli omejitve tipa double. To tudi pojasnjuje, zakaj kumulativne funkcije ne računate s seštevanjem dolgega repa gostotnih členov: vsak člen prinaša svojo napako zaokroževanja, napake se med potekom vrste kopičijo, regularizirana posebna funkcija pa se seštevku v celoti izogne tako, da namesto tega ovrednoti hitro konvergentno vrsto ali verižni ulomek.
Vrsta pod diagonalo, verižni ulomek nad njo
Rutina nepopolne gama funkcije sprejme eno odločitev, preden kar koli izračuna: primerja x z a + 1. Ta meja ni naključna. Razvoj v potenčno vrsto za P(a, x) hitro konvergira, ko je x majhen v primerjavi z a, in počasi (sčasoma neuporabno), ko je x velik. Verižni ulomek ima nasprotne lastnosti. Pogon zato uporablja potenčno vrsto za x pod a + 1 in Lentzov verižni ulomek za x enak ali nad a + 1, pri čemer se od vsake veje zahteva le delo, za katero je primerna.
Verižni ulomek potrebuje eno zaščito. Lentzova metoda deluje tako, da prenaša tekoči števec in imenovalec ter invertira imenovalec na vsakem koraku; če se kateri koli od njiju približa ničli, se inverzija zruši. Rešitev je majhna spodnja meja: kadarkoli vmesni člen pade pod velikost približno 1e-30, se omeji na 1e-30, kar ohranja rekurzijo končno, ne da bi motilo konvergirano vrednost. Enaka omejitev se iz istega razloga pojavi v verižnem ulomku nepopolne beta funkcije. Gre za majhno konstanto, ki opravlja pomembno delo – razliko med stabilnim vrednotenjem in deljenjem z nečim, kar je nerazločljivo od ničle.
Zgornji rep, Q(a, x), je preprosto 1 minus P(a, x) in tako se izračuna Poissonova kumulativna veja: verjetnost največ k dogodkov s srednjo vrednostjo λ je Q(k + 1, λ). Preusmeritev skozi zgornjo nepopolno gama funkcijo namesto seštevanja k + 1 Poissonovih členov je znova izbira za vrednotenje enega konvergentnega izraza namesto kopičenja številnih majhnih.
Diskretne mase brez prekoračitve fakultete
Diskretne porazdelitve prinašajo drugačno nevarnost. Binomska verjetnostna masa vključuje binomski koeficient, koeficient za "dvainpetdeset nad šestindvajset" pa je ogromno celo število. Če ga ustvarite neposredno, števec prekorači tip double še pred deljenjem, ki bi ga vrnilo na smiselno verjetnost. Pogon ga nikoli ne ustvari. Fakultete izračuna v logaritemskem prostoru prek log-gama funkcije, sešteje in odšteje logaritme, vključi logaritem verjetnosti uspeha in neuspeha ter na samem koncu enkrat potencira.
// 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));Sama log-gama funkcija je Lanczosev približek, ki je natančen na celotni pozitivni osi in poceni za vrednotenje. Ker se vsaka velika količina hrani kot logaritem do končnega klica Exp, je največje število, ki ga rutina sploh materializira, sama verjetnost, ki znaša največ ena. Poissonova funkcija mase sledi istemu receptu, pri čemer en sam log-gama člen predstavlja fakulteto v imenovalcu. Zaprte oblike so posebej obravnavane na robovih, kjer je p natanko nič ali ena, tako da koda nikoli ne pokliče Ln(0). HotXLS vrne 0.2460938 za BINOM.DIST(5,10,0.5,FALSE) in 0.6766764 for kumulativno POISSON.DIST(2,2,TRUE), kar se ujema z Excelom skozi vse izpisane številke.
Inverzi z zamejitvijo krivulje naprej
Inverzna funkcija porazdelitve sprašuje nasprotno: glede na verjetnost poišči vrednost, katere kumulativna porazdelitev (CDF) ji je enaka. Samo en inverz v tem naboru ima hitro neposredno formulo. NORM.S.INV, inverzna standardna normalna porazdelitev, uporablja Acklamov racionalni približek, par polinomskih razmerij, natančnih do približno natančnosti dvojne natančnosti (double) na celotnem območju, razdeljenih na osrednje območje in dva repa. Gre za vrednotenje zaprte oblike brez iteracije.
Drugi inverzi nimajo takšne formule, zato jih pogon obrača numerično. Odgovor zameji s spodnjo in zgornjo mejo, izbranima iz podpore porazdelitve, nato pa razpolavlja: ovrednoti CDF naprej na sredini, premakne tisto mejo, ki ohranja ciljno verjetnost zaprto, in ponavlja, dokler interval ni ozek. Za inverza gama in hi-kvadrat se zamejitev začne pri ničli in velikodušni zgornji oceni, zgrajeni iz oblike in merila, pri čemer podvoji zgornjo mejo, če verjetnost še ni zaokrožena. Inverz t zamejuje simetrične meje, ki se širijo navzven; inverz F pa se razpolavlja na ne-negativnem intervalu. Strošek je nekaj deset vrednotenj CDF na klic, kar je pri hitrosti preglednic neopazno, korist pa je v tem, da je vsak inverz natančno tako natančen kot funkcija naprej, ki jo obrača. Zato dvosmerna pot, kot je CHISQ.DIST(CHISQ.INV(0.7,5),5,TRUE), vrne 0.7 skoraj do potankosti.
Desetiški logaritem, ki se je skrival v repu
Tukaj je napaka, ki jo je vredno omeniti, saj je takšne vrste, ki dolgo preživi. Acklamova rutina za inverzno normalno porazdelitev ima tri veje. Široka osrednja veja, ki se uporablja vedno, ko se verjetnost nahaja med približno 0,025 in 0,975, požene vhod skozi polinomsko razmerje brez kakršnega koli logaritma. Obe veji na repih, za zelo majhne ali zelo velike verjetnosti, najprej vzameta logaritem vhoda, saj se rep obnaša kot kvadratni koren minus naravnega logaritma p.
Zgodnja različica veje za rep je vzela desetiški logaritem (logarithm base-10) tam, kjer bi moral biti naravni logaritem. Slednja se razlikujeta za konstantni faktor približno 2.30, zato so bili rezultati na repih napačni za skladen, precejšen odmik. Kljub temu je funkcija pri vsakem mimobežnem preverjanju delovala pravilno, saj se mimobežna preverjanja izvajajo na sredini. Klic NORM.S.INV(0.5) vrne nič, NORM.S.INV(0.975) vrne učbeniških 1.959964, oba pa tečeta skozi osrednji polinom, ki sploh ne kliče logaritma. Napaka se je pojavila šele, ko je verjetnost prešla v rep, na primer NORM.S.INV(0.001), ki mora vrniti -3.0902323, namesto tega pa je vrnil napačno vrednost zaradi razmerja med naravnim in desetiškim logaritmom. Vsaka funkcija, ki je v svojem repu odvisna od inverzne normale, vključno s pomočniki za intervale zaupanja, je podedovala enak odmik. Lekcija je vsakdanja in draga: funkcija z vejnato strukturo potrebuje testne točke znotraj vsake veje, saj bo pravilna pogosta pot zlahka prikrila pokvarjeno redko pot. Popravek je bil menjava enega samega znaka iz desetiškega logaritma v naravni logaritem, vrednosti na repih pa so se takoj poravnale z Excelovimi.
Predznak x določa rep t-porazdelitve
Studentova kumulativna t-funkcija vsebuje drobno posebnost, ki jo je enostavno obrniti narobe. Njena vrednost izhaja iz regularizirane nepopolne beta funkcije, ovrednotene pri df / (df + x²), toda ta vrednost beta predstavlja verjetnost v repu zunaj velikosti x in ne kumulativne verjetnosti do x. Simetrična oblika t-porazdelitve pomeni, da je pretvorba odvisna od tega, na kateri strani ničle leži x.
// 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 meanZa x nad ničlo je kumulativna verjetnost ena minus polovica simetričnega repa; za x pod ničlo je to polovica tega repa; pri ničli pa je natanko ena polovica. Če vrnete vrednost beta neposredno, poročate o napačni strani porazdelitve, kar pri katerem koli neničelnem x pomeni odstopanje za celotno telo krivulje. Različice z desnim repom in dvema repoma gradijo na isti veji, zato se klic T.DIST.2T(1,1) vrne kot 0.5, T.DIST(1,1,TRUE) pa kot 0.75, inverz T.INV pa se razpolavlja glede na ta popravljeni CDF, tako da se krog sklene.
Nič od tega ni vidno iz celice in to je pričakovani rezultat. Napišete formulo in preberete številko, ki se ujema z Excelom. Če pogon razširjate z lastno logiko, so mehanizmi registracije funkcije obravnavani v našem članku o razširitvi pogona formul s funkcijami po meri, način, kako formule dosegajo druge liste in poimenovana območja, pa je opisan v članku o definiranih imenih in formulah čez delovne liste. Vse to je na voljo znotraj komponente preglednic HotXLS spreadsheet component za Delphi in C++Builder, skupaj z vmesniki API za branje, pisanje, grafikone in oblikovanje, ki so obravnavani drugje na tem blogu.