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 tilnærme i hånden. En regnearkskomponent, der hævder at være kompatibel med Excel, skal gengive disse værdier til det sidste ciffer, som Excel viser, hvilket betyder, at de numeriske metoder skal gengives, ikke kun 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 kerne af specialfunktioner, grenbeslutningerne, der holder aritmetikken stabil, og én invers normalfordelingsfejl, der gemte sig i halen i lang tid, fordi det almindelige tilfælde aldrig berørte den fejlbehæftede linje
Ét kald i regnearket, halvtreds fordelinger bag det
Funktionerne spænder over de familier, som en statistikprojektmappe rækker ud efter. Der er normalfordelingsfamilien, NORM.DIST og NORM.S.DIST med deres inverse; gamma- og chi-i-anden-familien, 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; og inferenshjælperne som f.eks. CONFIDENCE.T og CONFIDENCE.NORM. Set fra kalderens plads er hver enkelt en simpel formel. Du indtaster input i celler, beder projektmappen om at evaluere og læser resultatet
var
wb: IXLSWorkbook;
sh: IXLSWorksheet;
begin
wb := TXLSWorkbook.Create;
try
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
finally
// wb is an interface, it is freed automatically when it goes out of scope,
// but we can release it explicitly if needed
end;
end;Calculate-metoden på projektmappen kompilerer og evaluerer en ad hoc-formel mod det aktive ark og leverer en Variant tilbage. En enkelt detalje spænder ben for folk i første forsøg: Formelparseren bag Calculate tager semikolon som sin argumentseparator, så det er =SUM(A1;B1), ikke =SUM(A1,B1). Gemte celleformler beholder Excels standard-komma. Den samme evaluator afsender enhver statistisk funktion nedenfor, så når først en af disse virker i Calculate, følger resten den 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 specialfunktioner: Den regulariserede nedre ufuldstændige gamma-funktion, skrevet P(a, x), og den regulariserede ufuldstændige beta-funktion, 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'en er gamma-CDF'en med formen df/2 og skalaen 2. Gamma-CDF'en is direkte P(a, x). De kumulative t-, F- og binomialfordelingsfunktioner er alle værdier af den regulariserede ufuldstændige beta-funktion ved de rigtige argumenter. Poissons CDF er den øvre ufuldstændige gamma Q. Implementer gamma- og betafunktionerne godt, og et dusin fordelinger arver deres nøjagtighed gratis
Ordet "regulariseret" er hele pointen. Den rå ufuldstændige gamma vokser som et fakultet, og det rå beta-integral kan underløbe eller overløbe længe før svaret gør det. De regulariserede former er divideret med den fuldstændige gamma eller beta, så de lever udelukkende i intervallet fra nul til et, hvilket er præcis det område, en sandsynlighed optager. Denne 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 man ikke beregner en CDF by adding up a long tail of density terms: hvert led bærer sin egen afrundingsfejl, fejl ophobes, efterhånden som rækken kører, og den regulariserede specialfunktion omgår summen helt ved i stedet at evaluere en hurtigt konvergerende række eller kædebrøk
Række 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. Denne grænse er ikke vilkårlig. Potensrækkens udvikling af P(a, x) konvergerer hurtigt, når x is lille i forhold til a, og langsomt, til sidst ubrugeligt, når x is stort. 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 gøre det arbejde, den er god til
Kædebrøken har brug for en enkelt beskyttelse. Lentz' metode fungerer ved at føre en løbende tæller og nævner og invertere nævneren på hvert trin, og hvis en af dem nærmer sig nul, sprænges inversionen i luften. Løsningen er en lige bundgrænse: Hver gang et mellemliggende led falder til under ca. 1e-30 i størrelse, fastklemmes det til 1e-30, hvilket holder rekursionen endelig uden at forstyrre den konvergerede værdi. Den samme fastklemning optræder i den ufuldstændige beta-funktions kædebrøk af samme årsag. Det er en lille konstant, der udfører et stort stykke arbejde, forskellen på 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 lede det gennem den øvre ufuldstændige gamma frem for at summere k + 1 Poisson-led er igen et valg om at evaluere ét konvergerende udtryk i stedet for at akkumulere mange små
Diskrete masser uden fakultetsoverløb
De diskrete fordelinger udgør en anden fare. En binomial sandsynlighedsmasse involverer en binomialkoefficient, og koefficienten for tooghalvtreds-vælg-seksogtyve er et enormt heltal. Hvis du danner det direkte, vil tælleren overløbe en double før den division, der ville bringe det tilbage til en fornuftig sandsynlighed. Motoren danner det aldrig direkte. Den beregner fakulteterne i log-rummet via log-gamma-funktionen, adderer og subtraherer logaritmerne, indfolder logaritmen for succes- og fejlsandsynlighederne og ophæver logaritmen (exponentierer) én enkelt 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 massefunktion følger den samme opskrift, med det enkelte log-gamma-led stående i stedet for fakultetet i nævneren. De lukkede former behandles som specialtilfælde ved grænserne, hvor p er præcis nul eller et, så koden kalder aldrig Ln(0). HotXLS returnerer 0.2460938 for BINOM.DIST(5,10,0.5,FALSE) og 0.6766764 for den kumulative POISSON.DIST(2,2,TRUE), hvilket matcher Excel i alle de viste cifre
Inverse funktioner 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 funktion i dette sæt har en hurtig direkte formel. NORM.S.INV, den inverse standardnormalfordeling, bruger en rationel Acklam-tilnærmelse, et par polynomielle forhold, der er nøjagtige til ca. præcisionen af en double over hele området, 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 svaret med en lav og høj grænse valgt ud fra fordelingens definitionsmængde, og foretager derefter bisektion: Evaluerer den fremadrettede CDF på midtpunktet, flytter den grænse, der holder målsandsynligheden indkredset, og gentager, 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 formen og skalaen, og fordobler den øvre grænse, hvis sandsynligheden endnu ikke er indkredset. Den inverse t-fordeling indkredser symmetriske grænser, der udvider sig udad; den inverse F-fordeling foretager bisektion på et ikke-negativt interval. Prisen er et par dusin CDF-evalueringer pr. kald, hvilket er usynligt ved regnearkshastighed, og fordelen er, at enhver invers funktion er præcis lige så nøjagtig som den fremadrettede funktion, den inverterer. Det er derfor, en tur-retur som f.eks. CHISQ.DIST(CHISQ.INV(0.7,5),5,TRUE) returnerer 0.7 på et hængende hår
Den base-10-logaritme, der gemte sig i halen
Her er fejlen, der er værd at fortælle om, fordi det er den slags, der overlever i lang tid. Acklams inverse normalfordelingsrutine 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 polynomielt forhold uden nogen logaritme overhovedet. De to halegrene, for meget små eller meget store sandsynligheder, tager hver 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 base-10-logaritme, hvor den naturlige logaritme hørte hjemme. De to adskiller sig med en konstant faktor på ca. 2.30, so the tail results were wrong by a consistent, sizable margin. Og alligevel så funktionen fin ud i ethvert overfladisk tjek, fordi overfladiske tjek 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 i nærheden af forholdet mellem naturlig logaritme og base-10. Enhver funktion, der afhænger af den inverse normalfordeling i sin hale, herunder konfidensinterval-hjælperne, arvede den samme skævhed. Læren er banal og dyr: En funktion med en forgrenet struktur har brug for testpunkter i hver 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 base-10-log til naturlig log, og haleværdierne faldt på plads svarende til Excels
Fortegnet for x bestemmer t-fordelingens hale
Students kumulative t-fordeling bærer på en deltalje, der er nem at få galt i halsen. Dens værdi kommer fra den regulariserede ufuldstændige beta evalueret ved df / (df + x²), men denne betavæ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å
var
ib: Double;
begin
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 mean
end;For 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 præcis en halv. Returner betaværdien direkte, og du rapporterer den forkerte side af fordelingen, forskudt med hele kurvens krop for ethvert x forskelligt fra nul. Varianterne for højre hale og to haler bygger på den samme gren, hvorfor T.DIST.2T(1,1) kommer tilbage som 0.5 og T.DIST(1,1,TRUE) som 0.75, og den inverse T.INV foretager bisektion mod denne korrigerede CDF, så turen frem og tilbage lukker
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 i at registrere 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 referencer på tværs af ark. Alt dette leveres i HotXLS spreadsheet component til Delphi og C++Builder, sammen med de API'er til læsning, skrivning, diagrammering og formatering, der er dækket andre steder på denne blog