Funcții statistice Excel în Delphi: NORM, CHISQ, BETA

Tastați =NORM.DIST(115,100,15,TRUE) într-o celulă și Excel returnează 0.8413447 fără nicio ceremonie. Apelul arată ca o simplă căutare. Nu este așa. În spatele acelui număr se află distribuția normală cumulativă, o integrală fără o formă închisă, iar în spatele CHISQ.INV.RT și BETA.DIST se află funcții speciale pe care o bibliotecă atentă trebuie să le evalueze, nu să le aproximeze manual. O componentă de foi de calcul care pretinde compatibilitate cu Excel trebuie să reproducă aceste valori până la ultima cifră afișată de Excel, ceea ce înseamnă reproducerea metodelor numerice, nu doar a numelor de funcții.

HotXLS implementează peste cincizeci de astfel de funcții statistice, iar munca care le face corecte este aproape complet invizibilă din bara de formule. Acesta este un tur al modului în care motorul le calculează: nucleul comun de funcții speciale, deciziile de ramificare care mențin aritmetica stabilă și o eroare de inversare normală care s-a ascuns în coada distribuției mult timp deoarece cazul comun nu a atins niciodată linia defectă.

Un singur apel de foaie de calcul, cincizeci de distribuții în spate

Funcțiile acoperă familiile pe care le folosește un registru de calcul statistic. Există familia normală, NORM.DIST și NORM.S.DIST cu inversele lor; familia gamma și chi-pătrat, GAMMA.DIST, CHISQ.DIST, CHISQ.DIST.RT, CHISQ.INV.RT; familia beta, BETA.DIST și BETA.INV; distribuțiile de eșantionare T.DIST, T.DIST.2T, F.DIST și F.INV; perechea discretă BINOM.DIST și POISSON.DIST; și asistenții de inferență cum ar fi CONFIDENCE.T și CONFIDENCE.NORM. Din perspectiva apelantului, fiecare este o singură formulă. Setați intrările în celule, cereți registrului de calcul să evalueze și citiți rezultatul.

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;

Metoda Calculate a registrului de calcul compilează și evaluează o formulă ad-hoc în raport cu foaia activă și returnează un Variant. Un detaliu îi încurcă pe utilizatori la prima încercare: parserul de formule din spatele Calculate folosește punct și virgulă ca separator de argumente, deci este =SUM(A1;B1), nu =SUM(A1,B1). Formulele de celule stocate păstrează virgula standard din Excel. Același evaluator apelează fiecare funcție statistică de mai jos, așa că odată ce una dintre acestea funcționează în Calculate, restul urmează aceeași cale.

Cele două funcții pe care se bazează totul

Majoritatea distribuțiilor cumulative din acest set nu sunt calculate prin însumarea sau integrarea propriilor definiții. Ele sunt calculate din două funcții speciale: funcția gamma incompletă inferioară regularizată, scrisă P(a, x), și funcția beta incompletă regularizată, scrisă I_x(a, b). Intern acestea sunt ajutoarele pe care se bazează dispatcherele, iar lanțul este scurt. CDF-ul chi-pătrat este CDF-ul gamma cu forma df/2 și scala 2. CDF-ul gamma este P(a, x) direct. Funcțiile cumulative t, F și binomială sunt toate valori ale funcției beta incomplete regularizate la argumentele corecte. CDF-ul Poisson este funcția gamma incompletă superioară Q. Implementați corect funcțiile gamma și beta și o duzină de distribuții vor moșteni acuratețea lor în mod gratuit.

Cuvântul "regularizată" este esența. Funcția gamma incompletă brută crește ca un factorial, iar integrala beta brută poate suferi depășire negativă (underflow) sau pozitivă (overflow) cu mult înainte ca răspunsul să o facă. Formele regularizate sunt împărțite la funcția gamma sau beta completă, astfel încât trăiesc în întregime în intervalul de la zero la unu, care este exact intervalul pe care îl ocupă o probabilitate. Acea normalizare este cea care permite aceleiași rutine să deservească un chi-pătrat cu două grade de libertate și unul cu două sute, fără ca termenii intermediari să depășească limitele unui double. De asemenea, explică de ce nu calculați un CDF prin adunarea unei cozi lungi de termeni de densitate: fiecare termen poartă propria eroare de rotunjire, erorile se acumulează pe măsură ce seria rulează, iar funcția specială regularizată evită complet suma prin evaluarea unei serii rapid convergente sau a unei fracții continue în schimb.

Serii sub diagonală, fracții continue deasupra ei

Rutina gamma incompletă ia o decizie înainte de a calcula ceva: compară x cu a + 1. Acea graniță nu este arbitrară. Expansiunea în serie de puteri a lui P(a, x) converge rapid când x este mic în raport cu a și lent, în cele din urmă inutil, când x este mare. Fracția continuă are caracterul opus. Astfel, motorul folosește seria de puteri pentru x sub a + 1 și o fracție continuă Lentz pentru x la sau peste a + 1, iar fiecare ramură este solicitată să facă doar munca la care este bună.

Fracția continuă are nevoie de o protecție. Metoda lui Lentz funcționează prin menținerea unui numărător și a unui numitor rulați și inversarea numitorului la fiecare pas, iar dacă oricare dintre aceștia se apropie de zero, inversarea eșuează. Soluția este o limită inferioară minusculă: ori de câte ori un termen intermediar scade sub aproximativ 1e-30 în magnitudine, este fixat la 1e-30, ceea ce menține recurența finită fără a perturba valoarea convergentă. Aceeași fixare aparece în fracția continuă a funcției beta incomplete din același motiv. Este o mică constantă care face o muncă importantă, diferența dintre o evaluare stabilă și o împărțire la ceva imposibil de distins de zero.

Coada superioară, Q(a, x), este pur și simplu 1 minus P(a, x), și așa se calculează ramura cumulativă a lui Poisson: probabilitatea a cel mult k evenimente cu media λ este Q(k + 1, λ). Rutarea acesteia prin funcția gamma incompletă superioară în loc de însumarea a k + 1 termeni Poisson este, din nou, o alegere de a evalua o singură expresie convergentă în loc de a acumula multe expresii mici.

Mase discrete fără depășire de factorial

Distribuțiile discrete ridică un pericol diferit. O masă de probabilitate binomială implică un coeficient binomial, iar coeficientul pentru cincizeci și doi luate câte douăzeci și șase este un întreg enorm. Generați-l direct și numărătorul va depăși un double înainte de împărțirea care l-ar aduce înapoi la o probabilitate rezonabilă. Motorul nu îl generează niciodată direct. Acesta calculează factorialele în spațiul logaritmic prin funcția log-gamma, adună și scade logaritmii, adăugă logaritmul probabilităților de succes și eșec și exponențiază o singură dată la sfârșit.

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

Funcția log-gamma în sine este o aproximare Lanczos, precisă pe întreaga axă pozitivă și ieftină de evaluat. Deoarece fiecare cantitate mare este păstrată ca logaritm al său până la Exp final, cel mai mare număr pe care rutina îl materializează vreodată este probabilitatea în sine, care este de cel mult unu. Funcția de masă a lui Poisson urmează aceeași rețetă, cu singurul termen log-gamma ținând locul factorialului în numitor. Formele închise sunt tratate ca cazuri speciale la margini, unde p este exact zero sau unu, astfel încât codul nu apelează niciodată Ln(0). HotXLS returnează 0.2460938 pentru BINOM.DIST(5,10,0.5,FALSE) și 0.6766764 pentru cumulativul POISSON.DIST(2,2,TRUE), potrivindu-se cu Excel prin cifrele pe care le afișează.

Inverse prin încadrarea curbei înainte

O funcție de distribuție inversă pune întrebarea opusă: dată fiind o probabilitate, găsiți valoarea al cărei CDF este egal cu ea. O singură inversă din acest set are o formulă directă rapidă. NORM.S.INV, inversa normalei standard, folosește o aproximare rațională Acklam, o pereche de rapoarte polinomiale precise până la aproximativ precizia unui double pe întregul interval, împărțită într-o regiune centrală și două cozi. Este o evaluare în formă închisă fără iterație.

Celelalte inverse nu au o astfel de formulă, so motorul le inversează numeric. Acesta încadrează (brackets) răspunsul cu o limită inferioară și una superioară alese din suportul distribuției, apoi caută prin bisecție: evaluează CDF-ul înainte la mijloc, deplasează limita care menține probabilitatea țintă inclusă și repetă până când intervalul este îngust. Pentru inversele gamma și chi-pătrat, intervalul începe de la zero și o estimare superioară generoasă construită din formă și scală, dublând limita superioară dacă probabilitatea nu este încă inclusă. Inversa t încadrează limite simetrice care se lărgesc spre exterior; inversa F caută prin bisecție pe un interval nenegativ. Costul este de câteva zeci de evaluări CDF per apel, ceea ce este invizibil la viteza foilor de calcul, iar avantajul este că fiecare inversă este la fel de precisă ca funcția înainte pe care o inversează. De aceea, o călătorie dus-întors cum ar fi CHISQ.DIST(CHISQ.INV(0.7,5),5,TRUE) returnează 0.7 cu o precizie milimetrică.

Logaritmul în baza 10 care s-a ascuns în coadă

Iată eroarea care merită povestită, deoarece este genul care supraviețuiește mult timp. Rutina de normală inversă Acklam are trei ramuri. Ramura centrală largă, utilizată ori de câte ori probabilitatea se situează între aproximativ 0.025 și 0.975, rulează intrarea printr-un raport polinomial fără niciun logaritm. Cele două ramuri de coadă, pentru probabilități foarte mici sau foarte mari, preiau fiecare mai întâi un logaritm al intrării, deoarece coada se comportă ca rădăcina pătrată din minus logaritmul natural al lui p.

O versiune timpurie a ramurii de coadă a preluat un logaritm în baza 10 acolo unde trebuia să fie logaritmul natural. Cele două diferă printr-un factor constant de aproximativ 2.30, așa că rezultatele pe coadă au fost greșește cu o marjă consistentă și considerabilă. Și totuși, funcția părea în regulă la fiecare verificare ocazională, deoarece verificările ocazionale au loc în mijloc. NORM.S.INV(0.5) este zero, NORM.S.INV(0.975) este clasicul 1.959964, iar ambele trec prin polinomul central care nu apelează deloc un logaritm. Eroarea a apărut doar când o probabilitate a trecut într-o coadă, cum ar fi NORM.S.INV(0.001), care trebuie să returneze -3.0902323 și în schimb a venit cu raportul dintre logaritmul natural și cel în baza 10 decalat. Orice funcție care depinde de normala inversă în coada sa, inclusiv ajutoarele pentru intervalul de încredere, a moștenit aceeași deviație. Lecția este banală și costisitoare: o funcție cu o structură de ramificare are nevoie de puncte de testare în fiecare ramură, deoarece o cale comună corectă va masca cu succes una rară defectă. Soluția a fost o schimbare a unui singur token de la logaritmul în baza 10 la logaritmul natural, iar valorile cozii s-au aliniat la cele ale Excel.

Semnul lui x decide coada distribuției t

Funcția cumulativă Student t poartă o subtilitate care este ușor de înțeles invers. Valoarea sa provine din funcția beta incompletă regularizată evaluată la df / (df + x²), dar acea valoare beta este probabilitatea din coada de dincolo de magnitudinea lui x, nu probabilitatea cumulativă până la x. Forma simetrică a distribuției t înseamnă că conversia depinde de care parte a lui zero se află x.

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

Pentru x mai mare ca zero, probabilitatea cumulativă este unu minus jumătate din coada simetrică; pentru x mai mic ca zero, este jumătate din acea coadă; la zero este exact o jumătate. Returnați valoarea beta direct și veți raporta partea greșită a distribuției, decalată cu întregul corp al curbei pentru orice x diferit de zero. Variantele de coadă dreaptă și cu două cozi se bazează pe aceeași ramură, motiv pentru care T.DIST.2T(1,1) se întoarce ca 0.5 și T.DIST(1,1,TRUE) ca 0.75, iar inversa T.INV face bisecție în raport cu acest CDF corectat, astfel încât călătoria dus-întors se închide.

Nimic din toate acestea nu este vizibil din celulă și acesta este rezultatul dorit. Scrieți o formulă și citiți un număr care este în concordanță cu Excel. Dacă extindeți motorul cu propria logică, mecanismul de înregistrare a unei funcții este acoperit în ghidul nostru despre motorul de formule și funcțiile personalizate, iar modul în care formulele accesează foi diferite și intervale numite este acoperit în articolul despre nume definite și formule între foi. Totul este inclus în componenta de foi de calcul HotXLS pentru Delphi și C++Builder, alături de API-urile de citire, scriere, diagrame și formatare acoperite în alte părți ale acestui blog.