Excel-statistikfunktioner i Delphi: NORM, CHISQ, BETA

Skriv =NORM.DIST(115,100,15,TRUE) i en cell och Excel returnerar 0.8413447 utan krångel. Anropet ser ut som en sökning. Det är det inte. Bakom det där talet döljer sig den kumulativa normalfördelningen, en integral utan stängd form, och bakom CHISQ.INV.RT och BETA.DIST sitter speciella funktioner som ett noggrant bibliotek måste utvärdera, inte uppskatta för hand. En kalkylbladskomponent som gör anspråk på Excel-kompatibilitet måste reproducera dessa värden till sista siffran Excel visar, vilket innebär att man måste reproducera 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, grenbesluten som håller aritmetiken stabil och en invers-normal bugg som gömde sig i svansen under lång tid eftersom det vanliga fallet aldrig rörde vid den trasiga raden

Ett kalkylbladsanrop, femtio fördelningar bakom det

Funktionerna spänner över de familjer en statistikarbetsbok sträcker sig efter. Det finns normalfördelningsfamiljen, NORM.DIST och NORM.S.DIST med deras inverser; gamma- och chi-två-familjen, GAMMA.DIST, CHISQ.DIST, CHISQ.DIST.RT, CHISQ.INV.RT; beta-familjen, BETA.DIST och BETA.INV; samplingsfördelningarna T.DIST, T.DIST.2T, F.DIST och F.INV; det diskreta paret BINOM.DIST och POISSON.DIST; samt inferenshjälpare som CONFIDENCE.T och CONFIDENCE.NORM. Från användarens sida är var och en av dem en enda formel. Du ställer in indata i celler, ber arbetsboken utvärdera och läser av 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;

Metoden Calculate på arbetsboken kompilerar och utvärderar en ad-hoc-formel mot det aktiva bladet och returnerar 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 skrivs =SUM(A1;B1), inte =SUM(A1,B1). Sparade cellformler behåller Excels standardkommatecken. Samma utvärderare skickar alla statistiska funktioner nedan, så så fort en av dessa fungerar i Calculate följer resten samma väg

De två funktionerna som allt annat bygger på

De flesta kumulativa fördelningar i den här uppsättningen 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 I_x(a, b). Internt är detta hjälparna som sändarna lutar sig mot, och kedjan är kort. Chi-två-fördelningens kumulativa fördelningsfunktion (CDF) är gamma-CDF:en med formen df/2 och skalan 2. Gamma-CDF:en är P(a, x) direkt. De kumulativa funktionerna för t, F och binomial är alla värden för den regulariserade ofullständiga beta-funktionen 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 helt gratis

Ordet "regulariserad" är hela poängen. Den råa ofullständiga gammafunktionen växer som en fakultet och den råa beta-integralen kan under- eller överflöda långt innan svaret gör det. De regulariserade formerna delas genom den fullständiga gamma- eller betafunktionen, så de lever helt i intervallet från noll till ett, vilket är exakt det intervall en sannolikhet upptar. Den normaliseringen är det som gör att samma rutin kan tjäna en chi-två med två frihetsgrader och en med tvåhundra utan att de mellanliggande termerna rinner iväg från slutet av 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 sin egen avrundningsfel, felen ackumuleras när serien körs, och den regulariserade specialfunktionen kringgår summan helt genom att istället utvärdera en snabbt konvergerande serie eller ett 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. Potensserie-expansionen av 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 vid eller över a + 1, och varje gren ombeds 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 litet golv: närhelst en mellanliggande term faller under cirka 1e-30 i storlek begränsas den till 1e-30, vilket håller rekursionen ändlig utan att störa det konvergerade värdet. Samma begränsning visas i den ofullständiga betafunktionens kedjebråk av samma anledning. Det är en liten konstant som gör ett tungt arbete, skillnaden mellan en stabil utvärdering och en division med något som inte kan skiljas från noll

Den övre svansen, Q(a, x), is 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 det via den övre ofullständiga gammafunktionen istället för att summera k + 1 Poisson-termer är, återigen, ett val att utvärdera ett konvergerande uttryck istället för att ackumulera many små

Diskreta massor utan spill av fakulteter

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 innan divisionen som skulle bringa den tillbaka till en vettig sannolikhet. Motorn skapar den aldrig. Den beräknar fakulteterna i log-rymd via log-gamma-funktionen, adderar och subtraherar logaritmerna, lägger till logaritmen för sannolikheten för framgång och misslyckande, och exponentierar en enda gång 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 kvantitet hålls som sin logaritmer fram till den sista Exp-anropet, är det största talet som rutinen någonsin skapar sannolikheten själv, vilken är högst ett. Poissons massfunktion följer samma recept, med den enskilda log-gamma-termen som ersätter fakulteten i nämnaren. De stängda formerna är specialbehandlade vid kanterna, där p är exakt noll eller ett, så koden anropar aldrig Ln(0). HotXLS returnerar 0.2460938 för BINOM.DIST(5,10,0.5,FALSE) och 0.6766764 för den kumulativa POISSON.DIST(2,2,TRUE), vilket matchar Excel genom alla siffror den skriver ut

Inverser genom att ringa in framåt-kurvan

En invers fördelningsfunktion ställer den motsatta frågan: 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, den inversa standardnormalfördelningen, använder en Acklam-rationell approximation, ett par polynomkvoter som är noggranna till ungefär en doubles precision över hela intervallet, uppdelad i en central region och två svansar. Det är en stängd utvärdering utan iteration

De andra inverserna har ingen sådan formel, så motorn inverterar dem numeriskt. Den ringar in svaret med en nedre och övre gräns vald från fördelningens stöd, och halverar sedan (bisection): 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 börjar inringningen på noll och en generös övre uppskattning byggd från form och skala, och fördubblar den övre gränsen om sannolikheten inte är innesluten ännu. T-inversen innesluter symmetriska gränser som vidgas utåt; F-inversen halverar på ett icke-negativt intervall. Kostnaden är ett dussintal CDF-utvärderingar per anrop, vilket är osynligt i kalkylbladshastighet, och fördelen är att varje invers är exakt lika noggrann som framåtfunktionen den inverterar. Det är därför en rundtur som CHISQ.DIST(CHISQ.INV(0.7,5),5,TRUE) returnerar 0.7 på ett hår när

Logaritmen med bas 10 som gömde sig i svansen

Här är den bugg som är värd att berätta, eftersom den är av den sort som överlever länge. Acklams invers-normala rutin har tre grenar. Den breda centrala grenen, som används närhelst sannolikheten ligger mellan cirka 0.025 och 0.975, kör indata genom en polynomkvot utan någon logaritmer någonstans i den. De två svansgrenarna, för mycket små eller mycket stora sannolikheter, tar båda en logaritmer av indata först, eftersom svansen beter sig som kvadratroten ur minus den naturliga logaritmen av p

En tidig version av svansgrenen tog en logaritmer med bas 10 där den naturliga logaritmen hörde hemma. De två skiljer sig åt med en konstant faktor på cirka 2.30, så svansresultaten blev fel med en konsekvent, betydande marginal. Ändå såg funktionen bra ut vid varje slumpmässig kontroll, eftersom enkla 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 och istället kom tillbaka nära den felaktiga kvoten för naturlig logaritmer kontra bas 10. Alla funktioner som beror på inversen av normalfördelningen i sin svans, inklusive hjälparna för konfidensintervall, ärvde samma skevhet. Läxan är vardaglig och dyr: en funktion med en grenstruktur behöver testpunkter inuti varje gren, eftersom en korrekt vanlig väg gladeligen maskerar en trasig ovanlig. Lösningen var en ändring av en enskild symbol från log10 till den naturliga logaritmen, och svansvärdena knäppte på plats till Excels

Tecknet på x avgör t-fördelningens svans

Student t-kumulativfunktionen bär på en finess som är lätt att få bakom bakfoten. Dess värde vemmer från den regulariserade ofullständiga betafunktionen utvärderad vid df / (df + x²), som det betavärdet är sannolikheten i svansen bortom storleken på x, inte den kumulativa sannolikheten upp till x. Den symmetriska formen på t-fördelningen innebär att omvandlingen beror på vilken sida om noll x hamnar

// 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 mean

Fö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 hälften. Returnera betavärdet direkt så rapporterar du fel sida av fördelningen, felaktigt med hela kurvans kropp för alla x som inte är noll. Varianten med högersvans och tvåsvans bygger på samma gren, vilket är anledningen till att T.DIST.2T(1,1) returnerar 0.5 och T.DIST(1,1,TRUE) returnerar 0.75, och den inversa T.INV halverar mot denna korrigerade CDF så att rundturen 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 beskrivs 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 beskrivs i artikeln om definierade namn och korsbladsreferenser. Allt levereras inuti HotXLS spreadsheet component för Delphi och C++Builder, tillsammans med API:er för läsning, skrivning, diagram och formatering som täcks på andra ställen i den här bloggen