Excel statistische functies in Delphi: NORM, CHISQ, BETA

Typ =NORM.DIST(115,100,15,TRUE) in een cel en Excel retourneert zonder poespas 0,8413447. De aanroep leest als een opvraging. Dat is het niet. Achter dat ene getal schuilt de cumulatieve normale verdeling, een integraal zonder gesloten vorm, en achter CHISQ.INV.RT en BETA.DIST schuilen speciale functies die een zorgvuldige bibliotheek moet evalueren, en niet met de hand benaderen. Een spreadsheetcomponent die compatibiliteit met Excel claimt, moet deze waarden reproduceren tot op het laatste cijfer dat Excel toont, wat betekent dat de numerieke methoden moeten worden gereproduceerd, en niet alleen de functienamen.

HotXLS implementeert meer dan vijftig van deze statistische functies, en het werk dat ze correct maakt is bijna volledig onzichtbaar vanaf de formulebalk. Dit is een rondleiding door hoe de engine ze berekent: de gedeelde kern van speciale functies, de splitsingsbeslissingen die de berekeningen stabiel houden, en één inverse-normale bug die zich lange tijd in de staart verborg omdat het algemene geval de beschadigde regel nooit raakte.

Eén werkbladaanroep, vijftig verdelingen erachter

De functies omvatten de families waar een statistisch werkboek naar grijpt. Er is de normale familie, NORM.DIST en NORM.S.DIST met hun inversen; the gamma- en chi-kwadraatfamilie, GAMMA.DIST, CHISQ.DIST, CHISQ.DIST.RT, CHISQ.INV.RT; de bètafamilie, BETA.DIST en BETA.INV; de steekproefverdelingen T.DIST, T.DIST.2T, F.DIST en F.INV; het discrete paar BINOM.DIST en POISSON.DIST; en de inferentiehelpers zoals CONFIDENCE.T en CONFIDENCE.NORM. Vanuit het perspectief van de aanroeper is elk een enkele formule. U stelt invoerwaarden in cellen in, vraagt het werkboek te evalueren en leest het resultaat.

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 Calculate-methode op het werkboek compileert en evalueert een ad-hocformule ten opzichte van het actieve blad en geeft een Variant terug. Eén detail zorgt bij de eerste poging vaak voor verwarring: de formuleparser achter Calculate gebruikt de puntkomma als scheidingsteken voor argumenten, dus het is =SUM(A1;B1), en niet =SUM(A1,B1). Opgeslagen celformules behouden de voor Excel standaard komma. Dezelfde evaluator handelt elke statistische functie hieronder af, dus zodra een van deze werkt in Calculate, de rest volgt hetzelfde pad.

De twee functies waarop al het andere is gebouwd

De meeste cumulatieve verdelingen in deze set worden niet berekend door hun eigen definities op te tellen of te integreren. Ze worden berekend op basis van twee speciale functies: de geregulariseerde lagere onvolledige gamma, geschreven als P(a, x), en de geregulariseerde onvolledige bèta, geschreven als I_x(a, b). Intern zijn dit de helpers waar de dispatchers op leunen, en de keten is kort. De chi-kwadraat-CDF is de gamma-CDF met vorm df/2 en schaal 2. De gamma-CDF is rechtstreeks P(a, x). De cumulatieve functies voor t, F en binomiaal zijn allemaal waarden van de geregulariseerde onvolledige bèta bij de juiste argumenten. De CDF van Poisson is de bovenste onvolledige gamma Q. Implementeer de gamma- en bètafuncties goed en een tiental verdelingen erft hun nauwkeurigheid gratis over.

Het woord 'geregulariseerd' is waar het om draait. De ruwe onvolledige gamma groeit als een faculteit en de ruwe bèta-integraal kan onderlopen (underflow) of overlopen (overflow) lang voordat het antwoord dat doet. De geregulariseerde vormen worden gedeeld door de volledige gamma of bèta, zodat ze volledig binnen het interval van nul tot één liggen, wat exact het bereik is dat een waarschijnlijkheid inneemt. Die normalisatie is wat ervoor zorgt dat dezelfde routine een chi-kwadraat met twee vrijheidsgraden kan bedienen en één met tweehonderd zonder dat de tussentijdse termen buiten het bereik van een double vallen. Het verklaart ook waarom u een CDF niet berekent door een lange staart van dichtheidstermen op te tellen: elke term brengt zijn eigen afrondingsfout met zich mee, de fouten stapelen zich op naarmate de reeks vordert, en de geregulariseerde speciale functie omzeilt de som volledig door in plaats daarvan een snel convergerende reeks of kettingbreuk te evalueren.

Reeks onder de diagonaal, kettingbreuk erboven

De onvolledige gamma-routine neemt één beslissing voordat ze iets berekent: ze vergelijkt x met a + 1. Die grens is niet willekeurig. De machtsreeksontwikkeling van P(a, x) convergeert snel wanneer x klein is ten opzichte van a, en langzaam, uiteindelijk nutteloos, wanneer x groot is. De kettingbreuk heeft het tegenovergestelde karakter. Dus de engine gebruikt de machtsreeks voor x onder a + 1 en een Lentz-kettingbreuk voor x op of boven a + 1, zodat elke tak alleen het werk hoeft te doen waar hij goed in is.

De kettingbreuk heeft één beveiliging nodig. De methode van Lentz werkt door een lopende teller en noemer bij te houden en de noemer bij elke stap te inverteren, en als een van beide nul naddert, de inversie mislukt. De oplossing is een minuscule ondergrens: telkens wanneer een tussentijdse term onder ongeveer 1e-30 in grootte zakt, wordt deze vastgeklemd (clamped) op 1e-30, wat de herhaling eindig houdt zonder de geconvergeerde waarde te verstoren. Dezelfde klem verschijnt om dezelfde reden in de kettingbreuk van de onvolledige bèta. Het is een kleine constante die belangrijk werk doet, het verschil tussen een stabile evaluatie en een deling door iets wat niet te onderscheiden is van nul.

De bovenstaart, Q(a, x), is simpelweg 1 min P(a, x), en dat is hoe de cumulatieve tak van Poisson wordt berekend: de waarschijnlijkheid van maximaal k gebeurtenissen met gemiddelde λ is Q(k + 1, λ). Het via de bovenste onvolledige gamma leiden in plaats van het optellen van k + 1 Poisson-termen is opnieuw een keuze om één convergerende expressie te evalueren in plaats van veel kleine te verzamelen.

Discrete massa's zonder faculteitsoverloop

De discrete verdelingen brengen een ander gevaar met zich mee. Een binomiale kansdichtheid (probability mass) omvat een binomiale coëfficiënt, en de coëfficiënt voor 'tweeënvijftig-kies-zesentwintig' is een enorm geheel getal. Bouw dit direct op en de teller loopt over (overflow) in een double nog vóór de deling die het terug zou brengen tot een zinnige waarschijnlijkheid. De engine bouwt dit nooit direct op. Hij berekent de faculteiten in de logaritmische ruimte via de log-gammafunctie, telt de logaritmen op en trekt ze af, voegt de logaritme van de succes- en faalkansen toe, en neemt pas helemaal aan het einde de e-macht.

// 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));

De log-gammafunctie zelf is een Lanczos-benadering, nauwkeurig over de hele positieve as en goedkoop te evalueren. Omdat elke grote hoeveelheid als logaritme wordt bewaard tot de uiteindelijke Exp, het grootste getal dat de routine ooit daadwerkelijk produceert is de waarschijnlijkheid zelf, die maximaal één is. De massafunctie van Poisson volgt hetzelfde recept, met de enkele log-gammaterm in de noemer de plaats inneemt van de faculteit. De gesloten vormen worden aan de randen als speciale gevallen behandeld, waar p exact nul of één is, zodat de code nooit Ln(0) aanroept. HotXLS retourneert 0,2460938 voor BINOM.DIST(5,10,0.5,FALSE) en 0,6766764 voor het cumulatieve POISSON.DIST(2,2,TRUE), wat overeenkomt met Excel tot in de getoonde cijfers.

Inversen door de voorwaartse curve in te sluiten

Een inverse verdelingsfunctie stelt de tegenovergestelde vraag: vind bij een gegeven waarschijnlijkheid de waarde waarvan de CDF eraan gelijk is. Slechts één inverse in deze set een snelle directe formule. NORM.S.INV, de inverse standaardnormale verdeling, maakt gebruik van een rationale benadering van Acklam, een paar polynoomverhoudingen die nauwkeurig zijn tot ruwweg de precisie van een double over het hele bereik, opgesplitst in een centraal gebied en twee staarten. Het is een evaluatie in gesloten vorm zonder iteratie.

De andere inversen hebben geen dergelijke formule, dus de engine inverteert ze numeriek. Hij sluit (bracket) het antwoord in met een onder- en bovengrens die zijn gekozen uit de drager van de verdeling, en halveert (bisecteert) vervolgens: evalueer de voorwaartse CDF op het middelpunt, verplaats de grens die de doelwaarschijnlijkheid ingesloten houdt, en herhaal dit tot het interval smal is. Voor de gamma- en chi-kwadraatinversen begint de insluiting bij nul en een genereuze bovengrens die is gebouwd op basis van de vorm en schaal, waarbij the bovengrens wordt verdubbeld als de waarschijnlijkheid nog niet is ingesloten. De t-inverse sluit symmetrische grenzen in die naar buiten toe breder worden; de F-inverse bisecteert op een niet-negatief interval. De kosten zijn een paar dozijn CDF-evaluaties per aanroep, wat onzichtbaar is op spreadsheetsnelheid, en het voordeel is dat elke inverse exact even nauwkeurig is als de voorwaartse functie die deze inverteert. Dat is waarom een heen- en terugreis zoals CHISQ.DIST(CHISQ.INV(0.7,5),5,TRUE) tot op een haar nauwkeurig 0,7 retourneert.

De logaritme met grondtal 10 die zich in de staart verbergt

Hier is de bug die het vertellen waard is, want het is het soort dat lang overleeft. De Acklam inverse-normale routine heeft drie takken. De brede centrale tak, die wordt gebruikt wanneer de waarschijnlijkheid tussen ongeveer 0,025 en 0,975 ligt, leidt de invoer door een polynoomverhouding zonder dat er ergens een logaritme in zit. De twee staarttakken, voor zeer kleine of zeer grote waarschijnlijkheden, nemen elk eerst een logaritme van de invoer, omdat de staart zich gedraagt als de vierkantswortel van minus de natuurlijke logaritme van p.

Een vroege versie van de staarttak nam een logaritme met grondtal 10 waar de natuurlijke logaritme hoorde. De twee verschillen met een constante factor van ongeveer 2,30, waardoor de resultaten in de staart met een consistente, aanzienlijke marge onjuist waren. En toch leek de functie prima in elke oppervlakkige controle, omdat oppervlakkige controles zich in het midden bevinden. NORM.S.INV(0.5) is nul, NORM.S.INV(0.975) is het klassieke 1,959964, en beide lopen door de centrale polynoom die helemaal geen logaritme aanroept. De error verscheen pas zodra een waarschijnlijkheid een staart binnenging, zeg NORM.S.INV(0.001), wat -3,0902323 moet retourneren en in plaats daarvan terugkwam met een afwijking door de verhouding tussen natuurlijke en grondtal-10 logaritmen. Elke functie die afhankelijk is van de inverse normale verdeling in zijn staart, inclusief de betrouwbaarheidsinterval-helpers, erfde dezelfde scheefheid. De les is alledaags en duur: een functie met een vertakkingsstructuur heeft testpunten nodig binnen elke tak, omdat een correct algemeen pad met plezier een kapot zeldzaam pad maskeert. De oplossing was een wijziging van één token van de grondtal-10 log naar de natuurlijke log, en de staartwaarden kwamen direct overeen met die van Excel.

Het teken van x bepaalt de staart van de t-verdeling

De Student t cumulatieve functie bevat een subtiliteit die gemakkelijk achterstevoren kan worden begrepen. De waarde ervan is afkomstig van de geregulariseerde onvolledige bèta geëvalueerd op df / (df + x²), maar die bètawaarde is de waarschijnlijkheid in de staart voorbij de grootte van x, en niet de cumulatieve waarschijnlijkheid tot aan x. De symmetrische vorm van de t-verdeling betekent dat de conversie afhangt van aan welke kant van nul x valt.

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

Voor x boven nul is de cumulatieve waarschijnlijkheid één min de helft van de symmetrische staart; voor x onder nul is het de helft van die staart; bij nul is het exact de helft. Retourneer de bètawaarde rechtstreeks en u rapporteert de verkeerde kant van de verdeling, met een afwijking van het volledige lichaam van de curve voor elke x die ongelijk is aan nul. De rechterstaart- en tweestaartvarianten bouwen voort op dezelfde tak, en dat is waarom T.DIST.2T(1,1) terugkomt als 0,5 en T.DIST(1,1,TRUE) als 0,75, en the inverse T.INV bisecteert ten opzichte van deze gecorrigeerde CDF zodat de cirkel rond is.

Niets hiervan is zichtbaar vanuit de cel, en dat is ook het beoogde resultaat. U schrijft een formule en leest een getal dat overeenkomt met Excel. Als u de engine uitbreidt met uw eigen logica, de mechanismen van het registreren van een functie worden behandeld in onze handleiding over de formule-engine en aangepaste functies. De manier waarop formules naar andere bladen en benoemde bereiken verwijzen, wordt behandeld in het artikel over gedefinieerde namen en formules over meerdere bladen. Dit alles wordt geleverd binnen de HotXLS-spreadsheetcomponent voor Delphi and C++Builder, naast de API's voor lezen, schrijven, grafieken en opmaak die elders op deze blog worden behandeld.