Skriv =NORM.DIST(115,100,15,TRUE) i en celle, og Excel returnerer 0,8413447 uden dikkedarer. Kaldet læses som et opslag. Det er det ikke. Bag det ene tal ligger den kumulative normalfordeling, et integral uden lukket form, og bag CHISQ.INV.RT og BETA.DIST sidder specielle funktioner, som et omhyggeligt bibliotek skal evaluere, ikke blot tilnærme i hånden. En regnearkskomponent, der hævder Excel-kompatibilitet, skal reproducere disse værdier til det sidste ciffer, Excel viser, hvilket betyder, at de numeriske metoder skal reproduceres, ikke blot funktionsnavnene.
HotXLS implementerer mere end halvtreds af disse statistiske funktioner, og det arbejde, der gør dem korrekte, er næsten helt usynligt fra formellinjen. Dette er en gennemgang af, hvordan motoren beregner dem: den fælles specielle funktionskerne, de forgreningsbeslutninger, der holder aritmetikken stabil, og en invers-normal-fejl, der gemte sig i halen i lang tid, fordi det almindelige tilfælde aldrig berørte den fejlbehæftede linje.
Et enkelt regnearkskald, halvtreds fordelinger bag det
Funktionerne spænder over de familier, som en statistikmappe rækker ud efter. Der er normalfordelingsfamilien, NORM.DIST og NORM.S.DIST med deres inverse; gamma- og chi-i-anden-familien (chi-square), GAMMA.DIST, CHISQ.DIST, CHISQ.DIST.RT, CHISQ.INV.RT; beta-familien, BETA.DIST og BETA.INV; stikprøvefordelingerne T.DIST, T.DIST.2T, F.DIST og F.INV; det diskrete par BINOM.DIST og POISSON.DIST; samt inferenshjælperne såsom CONFIDENCE.T og CONFIDENCE.NORM. Set fra kalderens plads er hver enkelt en enkelt formel. Du indsætter input i cellerne, beder projektmappen evaluere og læser resultatet.
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;Calculate-metoden på projektmappen kompilerer og evaluerer en ad-hoc formel mod det aktive ark og returnerer en Variant. Én detalje snyder folk i første forsøg: Formelparseren bag Calculate bruger semikolon som argumentseparator, så det er =SUM(A1;B1), ikke =SUM(A1,B1). Gemte celleformler beholder Excels standardkomma. Den samme evaluator afsender enhver statistisk funktion nedenfor, so når først en af disse fungerer i Calculate, følger resten samme sti.
De to funktioner, alt andet er bygget på
De fleste kumulative fordelinger i dette sæt beregnes ikke ved at summere eller integrere deres egne definitioner. De beregnes ud fra to specielle funktioner: den regulariserede nedre ufuldstændige gamma (regularized lower incomplete gamma), skrevet P(a, x), og den regulariserede ufuldstændige beta, skrevet Ix(a, b). Internt er disse de hjælpere, som afsenderne læner sig op ad, og kæden er kort. Chi-i-anden CDF er gamma-CDF med formen df/2 og skalaen 2. Gamma-CDF er P(a, x) direkte. De kumulative funktioner for t, F og binomial er alle værdier af den regulariserede ufuldstændige beta ved de rigtige argumenter. Poissons CDF er den øvre ufuldstændige gamma Q. Implementer gamma- og betafunktionerne godt, og en snes fordelinger arver deres nøjagtighed gratis.
Ordet "regulariseret" er hele pointen. Den rå ufuldstændige gamma vokser som en fakultet (factorial), og det rå beta-integral kan underløbe (underflow) eller overløbe (overflow) længe før svaret gør det. De regulariserede former er divideret med den fuldstændige gamma eller beta, så de lever helt i intervallet fra nul til et, hvilket er præcis det område, en sandsynlighed optager. Den normalisering er det, der lader den samme rutine betjene en chi-i-anden med to frihedsgrader og en med to hundrede, uden at de mellemliggende led løber tør for plads i en double. Det forklarer også, hvorfor du ikke beregner en CDF by adding up a long tail of density terms: Hvert led bærer sin egen afrundingsfejl, fejl akkumuleres, efterhånden som rækken kører, og den regulariserede specielle funktion omgår summen helt ved i stedet at evaluere en hurtigt konvergerende række eller et kædebrøk (continued fraction).
Rækker under diagonalen, kædebrøk over den
Den ufuldstændige gamma-rutine træffer én beslutning, før den beregner noget: Den sammenligner x med a + 1. Den grænse er ikke tilfældig. Potensrækkeudviklingen af P(a, x) konvergerer hurtigt, når x er lille i forhold til a, og langsomt – i sidste ende ubrugeligt – når x er stor. Kædebrøken har den modsatte karakter. Så motoren bruger potensrækken for x under a + 1 og en Lentz-kædebrøk for x på eller over a + 1, og hver gren bliver kun bedt om at udføre det arbejde, den er god til.
Kædebrøken har brug for en beskyttelse. Lentz' metode fungerer ved at føre en løbende tæller og nævner og inverterer nævneren ved hvert trin, og hvis en af dem nærmer sig nul, sprænges inversionen i luften. Løsningen er en lille bundgrænse (floor): Hver gang et mellemliggende led falder under ca. 1e-30 i størrelse, låses det fast på 1e-30, hvilket holder rekursionen endelig uden at forstyrre den konvergerede værdi. Den samme fastlåsning optræder i den ufuldstændige betas kædebrøk af samme årsag. Det er en lille konstant, der udfører et tungt arbejde – forskellen mellem en stabil evaluering og en division med noget, der ikke kan skelnes fra nul.
Den øvre hale, Q(a, x), er simpelthen 1 minus P(a, x), og det er sådan Poissons kumulative gren beregnes: Sandsynligheden for højst k hændelser med middelværdien λ er Q(k + 1, λ). At dirigere det gennem den øvre ufuldstændige gamma frem for at summere k + 1 Poisson-led er igen et valg om at evaluere ét konvergent udtryk i stedet for at akkumulere mange små.
Diskrete sandsynligheder uden fakultetsoverløb
De diskrete fordelinger medfører en anden fare. En binomial sandsynlighedsmasse involverer en binomialkoefficient, og koefficienten for tooghalvfjerds-vælg-seksogtyve (52 over 26) er et heltal. Hvis du danner det direkte, vil tælleren overløbe en double før den division, der ville bringe den tilbage til en fornuftig sandsynlighed. Motoren danner den aldrig direkte. Den beregner fakulteterne i log-rum via log-gamma-funktionen, lægger logaritmerne sammen og trækker dem fra hinanden, folder logaritmen for succes- og fiaskosandsynlighederne ind og eksponentierer én gang til allersidst.
// 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));Selve log-gamma-funktionen er en Lanczos-tilnærmelse, der er nøjagtig over hele den positive akse og billig at evaluere. Fordi enhver stor størrelse holdes som sin logaritme indtil den endelige Exp, er det største tal, rutinen nogensinde materialiserer, selve sandsynligheden, som højst er et. Poissons sandsynlighedsfunktion følger samme opskrift, hvor det enkelte log-gamma-led repræsenterer fakulteten i nævneren. De lukkede former behandles som specialtilfælde ved grænserne, hvor p er nøjagtig nul eller et, så koden aldrig kalder Ln(0). HotXLS returnerer 0,2460938 for BINOM.DIST(5,10,0.5,FALSE) and 0,6766764 for den kumulative POISSON.DIST(2,2,TRUE), hvilket matcher Excel i de cifre, den udskriver.
Inverse fordelinger ved at indkredse den fremadrettede kurve
En invers fordelingsfunktion stiller det modsatte spørgsmål: Givet en sandsynlighed, find den værdi, hvis CDF er lig med den. Kun én invers i dette sæt har en hurtig direkte formel. NORM.S.INV, den inverse standardnormalfordeling, bruger en rationel Acklam-tilnærmelse – et par polynomiske forhold, der er nøjagtige til omtrent præcisionen af en double over hele intervallet, opdelt i et centralt område og to haler. Det er en lukket form-evaluering uden iteration.
De andre inverse funktioner har ikke en sådan formel, så motoren inverterer dem numerisk. Den indkredser (bracket) svaret med en nedre og øvre grænse valgt fra fordelingens definitionsmængde og halverer derefter (bisect): Evaluer den fremadrettede CDF ved midtpunktet, flyt den grænse, der holder målsandsynligheden indkredset, og gentag, indtil intervallet er snævert. For gamma- og chi-i-anden-inverserne starter indkredsningen ved nul og et generøst øvre estimat bygget ud fra form og skala, idet den øvre grænse fordobles, hvis sandsynligheden endnu ikke er indkredset. Den inverse t indkredser symmetriske grænser, der udvider sig udad; den inverse F halverer på et ikke-negativt interval. Prisen er et par dusin CDF-evalueringer pr. kald, hvilket er usynligt ved regnearkshastighed, og fordelen er, at hver invers er præcis lige så nøjagtig som den fremadrettede funktion, den inverterer. Det er grunden til, at en rundturskonvertering som CHISQ.DIST(CHISQ.INV(0.7,5),5,TRUE) returnerer 0,7 inden for et hårstrå.
Den 10-talslogaritme, der gemte sig i halen
Her er fejlen, der er værd at fortælle, fordi det er den type, der overlever i lang tid. Acklam-invers-normal-rutinen har tre grene. Den brede centrale gren, der bruges, når sandsynligheden ligger mellem ca. 0,025 og 0,975, kører inputtet gennem et polynomisk forhold uden nogen logaritme overhovedet. De to halegrene, til meget små eller meget store sandsynligheder, tager hver især en logaritme af inputtet først, fordi halen opfører sig som kvadratroden af minus den naturlige logaritme af p.
En tidlig version af halegrenen tog en 10-talslogaritme, hvor den naturlige logaritme hørte hjemme. De to adskiller sig med en konstant faktor på ca. 2,30, så haleresultaterne var forkerte med en konsistent, betydelig margen. Og alligevel så funktionen fin ud i enhver tilfældig kontrol, fordi tilfældige kontroller lever i midten. NORM.S.INV(0.5) er nul, NORM.S.INV(0.975) er lærebogens 1,959964, og begge disse kører gennem det centrale polynomium, der slet ikke kalder en logaritme. Fejlen viste sig først, når en sandsynlighed krydsede ind i en hale, f.eks. NORM.S.INV(0.001), som skal returnere -3,0902323 og i stedet kom tilbage skæv med forholdet mellem den naturlige logaritme og 10-talslogaritmen. Enhver funktion, der afhænger af den inverse normalfordeling i sin hale, herunder hjælpefunktionerne til konfidensintervaller, arvede den samme skævhed. Læren er banal og dyr: En funktion med en forgreningsstruktur har brug for testpunkter inde i hver enkelt gren, fordi en korrekt fælles sti gladeligt vil maskere en ødelagt sjælden sti. Løsningen var en ændring af et enkelt symbol fra 10-tals log til den naturlige logaritme, og haleværdierne faldt på plads svarende til Excels.
Fortegnet for x bestemmer t-fordelingens hale
Student t-kumulative funktionen bærer en subtilitet, som er nem at få bagvendt. Dens værdi kommer fra den regulariserede ufuldstændige beta evalueret ved df / (df + x²), men den beta-værdi er sandsynligheden i halen ud over størrelsen af x, ikke den kumulative sandsynlighed op til x. Den symmetriske form af t-fordelingen betyder, at konverteringen afhænger af, hvilken side af nul x falder på.
// 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 meanFor x over nul er den kumulative sandsynlighed én minus halvdelen af den symmetriske hale; for x under nul er det halvdelen af denne hale; ved nul er det nøjagtig en halv. Returner beta-værdien direkte, og du rapporterer den forkerte side af fordelingen, forskudt med hele kurvens krop for ethvert x forskelligt fra nul. Højre-hale- og to-hale-varianterne bygger på samme gren, hvilket er grunden til, at T.DIST.2T(1,1) kommer tilbage som 0,5 og T.DIST(1,1,TRUE) som 0,75, og den inverse T.INV halverer mod denne korrigerede CDF, så rundturen lukkes.
Intet af dette er synligt fra cellen, og det er det tilsigtede resultat. Du skriver en formel og læser et tal, der stemmer overens med Excel. Hvis du udvider motoren med din egen logik, er mekanikken til registrering af en funktion dækket i vores gennemgang af formelmotoren og brugerdefinerede funktioner, og måden formler rækker på tværs af ark og navngivne områder er dækket i artiklen om definerede navne og formler på tværs af ark. Alt dette leveres inde i HotXLS regnearkskomponent til Delphi og C++Builder sammen med API'erne til læsning, skrivning, diagrammer og formatering, der er dækket andre steder på denne blog.