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 zoekactie. Dat is het niet. Achter dat ene getal gaat de cumulatieve normale verdeling schuil, een integraal zonder gesloten vorm, en achter CHISQ.INV.RT en BETA.DIST liggen speciale functies dat een zorgvuldige bibliotheek moet evalueren en niet met de hand mag benaderen. Een spreadsheetcomponent die claimt compatibel te zijn met Excel, 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 de manier waarop de engine ze berekent: de gedeelde kern van speciale functies, de vertakkingsbeslissingen die de rekenkunde stabiel houden, en één inverse-normaal-bug die lange tijd in de staart verborgen bleef omdat het veelvoorkomende geval de kapotte regel nooit raakte

Eén werkbladaanroep, vijftig verdelingen erachter

De functies omvatten de families waarnaar een statistisch werkboek grijpt. Er is de normale familie, NORM.DIST en NORM.S.DIST met hun inversen; de gamma- en chi-kwadraat-familie, GAMMA.DIST, CHISQ.DIST, CHISQ.DIST.RT, CHISQ.INV.RT; de bèta-familie, 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. Vanaf de stoel van de aanroeper is elk daarvan een enkele formule. U voert invoergegevens in cellen in, vraagt de werkmap om deze te evalueren en leest het resultaat

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.Free is handled by interfaces
  end;
end;

De methode Calculate op de werkmap compileert en evalueert een ad-hocformule tegen het live-blad en geeft een Variant terug. Eén detail brengt mensen bij de eerste poging in verwarring: de formule-parser achter Calculate gebruikt de puntkomma als argument-scheidingsteken, dus het is =SUM(A1;B1), niet =SUM(A1,B1). Opgeslagen celformules behouden de Excel-standaard komma. Dezelfde evaluator handelt elke hieronder beschreven statistische functie af, dus zodra een daarvan werkt in Calculate, volgt de rest hetzelfde pad

De twee functie 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 onderste onvolledige gamma, geschreven als P(a, x), en de geregulariseerde onvolledige bèta, geschreven als I_x(a, b). Intern zijn dit de hulpprogramma's 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 direct P(a, x). De t-, F- en binomiale cumulatieve functies zijn allemaal waarden van de geregulariseerde onvolledige bèta bij de juiste argumenten. Poisson's CDF is de bovenste onvolledige gamma Q. Implementeer de gamma- en bèta-functies goed en een dozijn verdelingen erft hun nauwkeurigheid gratis over

Het woord "geregulariseerd" (regularized) is waar het om draait. De ruwe onvolledige gamma groeit als een faculteit en de ruwe bèta-integraal kan underflowen of overflowen lang voordat het antwoord dat doet. De geregulariseerde vormen worden gedeeld door de volledige gamma of bèta, zodat ze volledig in het interval van nul tot één liggen, wat precies 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 een met tweehonderd, zonder dat de tussentijdse termen over de limiet van een double lopen. Het verklaart ook waarom u een CDF niet berekent door een lange staart van dichtheidstermen op te tellen: elke term draagt 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

Reeksen onder de diagonaal, kettingbreuk erboven

De onvolledige gamma-routine neemt één beslissing voordat ze berekent: ze vergelijkt x met a + 1. Die grens is niet willekeurig. De machtsreeksontwikkeling van P(a, x) converteert 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, en elke tak wordt gevraagd alleen het werk te doen waar ze goed in is

De kettingbreuk heeft één beveiliging nodig. De methode van Lentz werkt door een lopende teller en noemer mee te voeren en de noemer bij elke stap om te keren, en als een van beide de nul naderen, explodeert de omkering. De oplossing is een minuscule ondergrens: telkens wanneer een tussentijdse term onder ongeveer 1e-30 in grootte zakt, wordt deze vastgezet 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 zwaar werk verricht: het verschil tussen een stabiele evaluatie en een deling door iets dat niet van nul te onderscheiden is

De bovenste staart, Q(a, x), is simpelweg 1 minus P(a, x), en dat is hoe de cumulatieve tak van Poisson wordt berekend: de kans op maximaal k gebeurtenissen met gemiddelde λ is Q(k + 1, λ). Het leiden hiervan via de bovenste onvolledige gamma in plaats van het optellen van k + 1 Poisson-termen is, nogmaals, een keuze om één convergerende uitdrukking te evalueren in plaats van veel kleine uitdrukkingen te accumuleren

Discrete massa's zonder faculteit-overflow

De discrete verdelingen brengen een ander gevaar met zich mee. Een binomiale kansmassa omvat een binomiale coëfficiënt, en de coëfficiënt voor 52-boven-26 is een enorm geheel getal. Vorm het rechtstreeks en de teller overschrijdt de limiet van een double nog voor de deling die het terug zou brengen tot een zinvolle kans. De engine vormt dit nooit rechtstreeks. Ze berekent de faculteiten in log-ruimte via de log-gamma-functie, telt de logs op en trekt ze af, voegt de log van de succes- en faalkansen toe, en neemt pas helemaal aan het einde de exponent

// 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-gamma-functie zelf is een Lanczos-benadering, nauwkeurig over de hele positieve as en goedkoop te evalueren. Omdat elke grote hoeveelheid als zijn logaritme wordt vastgehouden tot de uiteindelijke Exp, is het grootste getal dat de routine ooit materialiseert de kans zelf, die maximaal één is. Poisson's massafunctie volgt hetzelfde recept, waarbij de enkele log-gamma-term staat voor de faculteit in de noemer. De gesloten vormen zijn speciaal behandeld aan de randen, 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 de cumulatieve POISSON.DIST(2,2,TRUE), wat overeenkomt met Excel tot op de cijfers die het afdrukt

Inversen door de voorwaartse curve in te sluiten (bracketen)

Een inverse verdelingsfunctie stelt de tegenovergestelde vraag: vind bij een gegeven kans de waarde waarvan de CDF eraan gelijk is. Slechts één inverse in deze set heeft een snelle, directe formule. NORM.S.INV, de inverse standaardnormaal, gebruikt een rationele benadering van Acklam, een paar polynoomverhoudingen die tot ongeveer de precisie van een double nauwkeurig zijn over het hele bereik, verdeeld 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. Ze sluit het antwoord in met een onder- en bovengrens gekozen uit de ondersteuning van de verdeling, en halveert vervolgens (bisection): evalueer de voorwaartse CDF op het middelpunt, verplaats de grens die de doelkans ingesloten houdt, en herhaal dit totdat het interval smal is. Voor de gamma- en chi-kwadraatinversen begint het insluiten bij nul en een genereuze bovengrens schatting gebouwd uit de vorm en schaal, waarbij de bovengrens wordt verdubbeld als de kans nog niet is ingesloten. De t-inverse sluit symmetrische grenzen in die naar buiten toe breder worden; de F-inverse halveert 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 ze inverteert. Dat is de reden waarom een round-trip zoals CHISQ.DIST(CHISQ.INV(0.7,5),5,TRUE) 0,7 retourneert tot op een haar nauwkeurig

De base-10-logaritme die in de staart verborgen zat

De Acklam inverse-normaal-routine heeft drie takken. De brede centrale tak, die wordt gebruikt wanneer de kans tussen ongeveer 0,025 en 0,975 ligt, leidt de invoer door een polynoomverhouding zonder dat er ergens een logaritme in voorkomt. De twee staarttakken, voor zeer kleine of zeer grote kansen, 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 base-10-logaritme waar de natuurlijke logaritme hoorde te staan. De twee verschillen met een constante factor van ongeveer 2,30, waardoor de resultaten in de staart consistent en aanzienlijk afweken. En toch leek de functie prima bij elke vluchtige controle, omdat vluchtige controles zich in het midden bevinden. NORM.S.INV(0.5) is nul, NORM.S.INV(0.975) is de bekende 1,959964, en beide lopen door de centrale polynoom die helemaal geen logaritme aanroept. De fout verscheen pas zodra een kans in een staart terechtkwam, bijvoorbeeld NORM.S.INV(0.001), wat -3,0902323 moet retourneren en in plaats daarvan terugkwam met een afwijking die dicht bij de verhouding tussen de natuurlijke en base-10-logaritme lag. Elke functie die in zijn staart afhankelijk is van de inverse standaardnormaal, inclusief de betrouwbaarheidsintervalhelpers, erfde dezelfde afwijking. De les is alledaags en duur: een functie met een vertakkingsstructuur heeft testpunten nodig binnen elke tak, omdat een correct algemeen pad vrolijk een kapot zeldzaam pad zal maskeren. De oplossing was een verandering van één token van de base-10-log naar de natuurlijke log, en de staartwaarden kwamen direct overeen met die van Excel

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

De cumulatieve Student-t-verdeling bevat een subtiliteit die men makkelijk verkeerd begrijpt. De waarde ervan is afkomstig van de geregulariseerde onvolledige bèta geëvalueerd op df / (df + x²), maar die bèta-waarde is de kans in de staart voorbij de grootte van x, niet de cumulatieve kans 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 zero is de cumulatieve kans één minus de helft van de symmetrische staart; voor x onder zero is dat de helft van die staart; bij zero is het exact de helft. Retourneer de bèta-waarde direct en u rapporteert de verkeerde kant van de verdeling, die met het hele lichaam van de curve afwijkt voor elke x die ongelijk is aan zero. De rechtszijdige en tweezijdige varianten bouwen voort op dezelfde tak, wat verklaart waarom T.DIST.2T(1,1) terugkomt als 0,5 en T.DIST(1,1,TRUE) als 0,75, en de inverse T.INV halveert tegen deze gecorrigeerde CDF zodat de round-trip klopt

Niets van dit alles 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, worden de mechanismen voor het registreren van een functie behandeld in onze handleiding over de formule-engine en aangepaste functies, en de manier waarop formules over bladen en benoemde bereiken heen reiken, wordt behandeld in het artikel over gedefinieerde namen en cross-sheet-verwijzingen. Dit alles wordt geleverd in de HotXLS-spreadsheetcomponent voor Delphi en C++Builder, naast de API's voor lezen, schrijven, grafieken maken en formatteren die elders op dit blog worden behandeld