Skriv =NORM.DIST(115,100,15,TRUE) i en cell och Excel returnerar 0.8413447 utan omsvep. Anropet kan läsas som en uppslagning. Det är det inte. Bakom det där talet finns den kumulativa normalfördelningen, en integral utan sluten form, och bakom CHISQ.INV.RT och BETA.DIST sitter speciella funktioner som ett noggrant bibliotek måste utvärdera, inte bara uppskatta för hand. En kalkylblads-komponent som hävdar Excel-kompatibilitet måste återskapa dessa värden till sista siffran Excel visar, vilket innebär att återskapa de numeriska metoderna, inte bara funktionsnamnen.
HotXLS implementerar mer än femtio av dessa statistiska funktioner, och arbetet som gör dem korrekta är nästan helt osynligt från formelfältet. Detta är en rundtur i hur motorn beräknar dem: den delade kärnan för specialfunktioner, de förgreningsbeslut som håller aritmetiken stabil, och en invers normalfördelningsbugg som gömde sig i svansen under en lång tid eftersom det vanliga fallet aldrig rörde den trasiga kodraden.
Ett kalkylbladsanrop, femtio fördelningar bakom det
Funktionerna spänner över de familjer en statistisk arbetsbok sträcker sig efter. Det finns normalfördelningsfamiljen, NORM.DIST och NORM.S.DIST med sina inverser; gamma- och chi-två-familjen, GAMMA.DIST, CHISQ.DIST, CHISQ.DIST.RT, CHISQ.INV.RT; betafamiljen, BETA.DIST och BETA.INV; urvalsfördelningarna T.DIST, T.DIST.2T, F.DIST och F.INV; de diskreta paren BINOM.DIST och POISSON.DIST; samt inferenshjälparna som CONFIDENCE.T och CONFIDENCE.NORM. Från anroparens sida är varje en enskild formel. Du ställer in indata i celler, ber arbetsboken utvärdera och läser resultatet.
Metoden Calculate på arbetsboken kompilerar och utvärderar en ad-hoc-formel mot det aktiva bladet och lämnar tillbaka en Variant. En detalj ställer till det för folk vid första försöket: formelparsern bakom Calculate använder semikolon som argumentseparator, så det är =SUM(A1;B1), inte =SUM(A1,B1). Sparade cellformler behåller Excels standardkommatecken. Samma utvärderare skickar vidare varje statistisk funktion nedan, så när en av dessa väl fungerar i Calculate följer resten samma väg.
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;De två funktioner som allt annat bygger på
De flesta kumulativa fördelningar i denna uppsättning beräknas inte genom att summera eller integrera sina egna definitioner. De beräknas från två specialfunktioner: den regulariserade nedre ofullständiga gammafunktionen, skriven P(a, x), och den regulariserade ofullständiga betafunktionen, skriven Ix(a, b). Internt är det dessa hjälpare som fördelarna lutar sig mot, och kedjan är kort. Chi-två-fördelningens kumulativa fördelningsfunktion (CDF) är gammafördelningens CDF med formen df/2 och skalan 2. Gammafördelningens CDF är P(a, x) direkt. De kumulativa funktionerna för t, F och binomialfördelningarna är alla värden av den regulariserade ofullständiga betafunktionen vid rätt argument. Poissons CDF är den övre ofullständiga gammafunktionen Q. Implementera gamma- och betafunktionerna väl så ärver ett dussin fördelningar deras noggrannhet gratis.
Ordet "regulariserad" är hela poängen. Den råa ofullständiga gammafunktionen växer som en fakultet och den råa betaintegralen kan understiga (underflow) eller överstiga (overflow) långt innan svaret gör det. De regulariserade formerna delas med den gamma- eller betafunktionen, så de lever helt i intervallet från noll till ett, vilket är exakt det område en sannolikhet upptar. Den normaliseringen är det som gör att samma rutin kan tjäna en chi-två-fördelning med två frihetsgrader och en med tvåhundra utan att de mellanliggande termerna rinner utanför gränsen för en double. Det förklarar också varför du inte beräknar en CDF genom att addera en lång svans av täthets-termer: varje term bär på sitt eget avrundningsfel, felen ackumuleras allt eftersom serien körs, och den regulariserade specialfunktionen kringgår summan helt genom att istället utvärdera en snabbt konvergerande serie eller kedjebråk.
Serier under diagonalen, kedjebråk över den
Rutinen för ofullständig gamma fattar ett beslut innan den beräknar något: den jämför x med a + 1. Den gränsen är inte godtycklig. Potensserieutvecklingen för P(a, x) konvergerar snabbt när x är litet i förhållande till a, och långsamt, slutligen oanvändbart, när x är stort. Kedjebråket har motsatt karaktär. Så motorn använder potensserien för x under a + 1 och ett Lentz-kedjebråk för x på eller över a + 1, och varje gren blir därmed ombedd att göra endast det arbete den är bra på.
Kedjebråket behöver ett skydd. Lentz metod fungerar genom att bära en löpande täljare och nämnare och invertera nämnaren vid varje steg, och om någon av dem närmar sig noll exploderar inversionen. Lösningen är ett pyttelitet golv: närhelst en mellanliggande term faller under ungefär 1e-30 i storlek begränsas den till 1e-30, vilket håller rekurrensen ändlig utan att störa det konvergerade värdet. Samma begränsning dyker upp i den ofullständiga betafunktionens kedjebråk av samma anledning. Det är en liten konstant som gör ett bärande arbete, skillnaden mellan en stabil utvärdering och en division med något som är omöjligt att skilja från noll.
Den övre svansen, Q(a, x), är helt enkelt 1 minus P(a, x), och det är så Poissons kumulativa gren beräknas: sannolikheten för högst k händelser med medelvärdet λ är Q(k + 1, λ). Att dirigera den genom den övre ofullständiga gammafunktionen snarare än att summera k + 1 Poisson-termer är, återigen, ett val att utvärdera ett konvergerande uttryck istället för att ackumulera många små.
Diskreta massor utan fakultetsspill
De diskreta fördelningarna medför en annan risk. En binomial sannolikhetsmassa involverar en binomialkoefficient, och koefficienten för femtiotvå-välj-tjugosex är ett enormt heltal. Skapa den direkt och täljaren spiller över en double före divisionen som skulle föra den tillbaka till en rimlig sannolikhet. Motorn skapar den aldrig. Den beräknar fakulteterna i log-rymden via log-gamma-funktionen, adderar och subtraherar logaritmerna, lägger till logaritmen för framgångs- och misslyckandesannolikheterna, och tar exponenten en gång helt i slutet.
// 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));Själva log-gamma-funktionen är en Lanczos-approximation, noggrann över hela den positiva axeln och billig att utvärdera. Eftersom varje stor storhet hålls som sin logaritmer fram till den sista Exp-funktionen, är det största talet rutinen någonsin materialiserar själva sannolikheten, vilken är högst ett. Poissons massfunktion följer samma recept, där den enskilda log-gamma-termen står för fakulteten i nämnaren. De slutna formerna behandlas som specialfall vid gränserna, där p är exakt noll eller ett, så att koden never calls Ln(0). HotXLS returnerar 0.2460938 för BINOM.DIST(5,10,0.5,FALSE) och 0.6766764 för the cumulative POISSON.DIST(2,2,TRUE), matching Excel through the digits it prints.
Inverser genom att ringa in framåt-kurvan
En invers fördelningsfunktion ställer motsatt fråga: givet en sannolikhet, hitta värdet vars CDF är lika med den. Endast en invers i denna uppsättning har en snabb direkt formel. NORM.S.INV, inversen av standardnormalfördelningen, använder en Acklam-rationell approximation, ett par polynomkvoter som är noggranna till ungefär precisionen av en double över hela intervallet, uppdelad i en central region och två svansar. Det är en utvärdering i sluten form utan iteration.
De andra inverserna har ingen sådan formel, så motorn inverterar dem numeriskt. Den ringar in (bracket) svaret med en nedre och övre gräns vald från fördelningens definitionsmängd (support), och halverar sedan (bisects): utvärdera framåt-CDF vid mittpunkten, flytta den gräns som håller målsannolikheten innesluten, och upprepa tills intervallet är smalt. För gamma- och chi-två-inverserna startar inringningen vid noll och en generös övre uppskattning byggd på form och skala, och fördubblar den övre gränsen om sannolikheten ännu inte är innesluten. T-inversen inringar symmetriska gränser som vidgas utåt; F-inversen halverar på ett icke-negativt intervall. Kostnaden är några dussin CDF-utvärderingar per anrop, vilket är osynligt i kalkylbladshastighet, och fördelen är att varje invers är exakt lika noggrann som den framåtfunktion den inverterar. Det är därför en rundresa som CHISQ.DIST(CHISQ.INV(0.7,5),5,TRUE) returnerar 0.7 så när som på ett hårstrå.
Den 10-baserade logaritmen som gömde sig i svansen
Här är buggen som är värd att berätta, eftersom den är av den typen som överlever under lång tid. Acklams rutin för invers normalfördelning har tre grenar. Den breda centrala grenen, som används närhelst sannolikheten ligger mellan ungefär 0.025 och 0.975, kör indatan genom en polynomkvot utan någon logaritmer i sig. De två svansgrenarna, för mycket små eller mycket stora sannolikheter, tar vardera en logaritmer av indatan först, eftersom svansen beter sig som kvadratroten av minus den naturliga logaritmen för p.
En tidig version av svansgrenen tog en 10-baserad logaritmer där den naturliga logaritmen hörde hemma. De två skiljer sig åt med en konstant faktor på ungefär 2.30, så svansresultaten var felaktiga med en konsekvent, betydande marginal. Och ändå såg funktionen bra ut i alla vardagliga kontroller, eftersom vardagliga kontroller lever i mitten. NORM.S.INV(0.5) är noll, NORM.S.INV(0.975) är lärobokens 1.959964, och båda dessa körs genom det centrala polynomet som aldrig anropar en logaritmer överhuvudtaget. Felet visade sig först när en sannolikhet korsade in i en svans, till exempel NORM.S.INV(0.001), som måste returnera -3.0902323 men istället kom tillbaka felaktig på grund av förhållandet mellan den naturliga och 10-baserade logaritmen. Alla funktioner som beror på den inversa normalfördelningen i sin svans, inklusive hjälparna för konfidensintervall, ärvde samma skevhet. Lärdomen är vardaglig och dyrbar: en funktion med en förgrenad struktur behöver testpunkter i varje gren, eftersom en korrekt vanlig sökväg glatt maskerar en trasig sällsynt. Lösningen var en förändring av en enda token från 10-baserad logg till naturlig logg, och svansvärdena smälte samman med Excels.
Tecknet för x avgör t-fördelningens svans
Student t-fördelningens kumulativa funktion bär på en subtilitet som är lätt att få bakom bakfoten. Dess värde kommer från den regulariserade ofullständiga betafunktionen utvärderad vid df / (df + x²), men det betavärdet är sannolikheten i svansen bortom storleken av x, inte den kumulativa sannolikheten upp till x. Den symmetriska formen av t-fördelningen innebär att konverteringen beror på vilken sida av noll som x faller.
// 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 meanFör x över noll är den kumulativa sannolikheten ett minus hälften av den symmetriska svansen; för x under noll är den hälften av den svansen; vid noll är den exakt en halv. Returnerar du betavärdet direkt rapporterar du fel sida av fördelningen, felaktig med hela kurvans kropp för varje nollskilt x. Varianterna för högersvans och två svansar bygger på samma gren, vilket är anledningen till att T.DIST.2T(1,1) kommer tillbaka som 0.5 och T.DIST(1,1,TRUE) som 0.75, och inversen T.INV halverar mot denna korrigerade CDF så att rundresan stängs.
Inget av detta är synligt från cellen, och det är det avsedda resultatet. Du skriver en formel och läser ett tal som stämmer överens med Excel. Om du utökar motorn med din egen logik täcks mekaniken för att registrera en funktion i vår genomgång av formelmotorn och anpassade funktioner, och hur formler når över blad och namngivna områden täcks i artikeln om definierade namn och formler över flera blad. Allt levereras inuti HotXLS spreadsheet component för Delphi och C++Builder, tillsammans med läs-, skriv-, diagram- och formaterings-API:erna som beskrivs på andra ställen i den här bloggen.