Funzioni statistiche di Excel in Delphi: NORM, CHISQ, BETA

Inserendo =NORM.DIST(115,100,15,TRUE) in una cella, Excel restituisce 0.8413447. La chiamata appare come una semplice ricerca. Ma non è così. Dietro quel singolo numero risiede la distribuzione normale cumulativa, un integrale privo di forma chiusa, e dietro CHISQ.INV.RT e BETA.DIST si trovano funzioni speciali che una libreria accurata deve calcolare effettivamente, senza ricorrere ad approssimazioni manuali. Un componente per fogli di calcolo che dichiara la compatibilità con Excel deve riprodurre questi valori fino all'ultima cifra mostrata da Excel, il che implica riprodurre i metodi numerici e non solo i nomi delle funzioni.

HotXLS implementa più di cinquanta di queste funzioni statistiche, e il lavoro che ne garantisce la correttezza è quasi del tutto invisibile dalla barra delle formule. Questo è un percorso di analisi del modo in cui il motore le calcola: il nucleo comune di funzioni speciali, le scelte di ramificazione che mantengono stabile il calcolo e un bug della normale inversa che è rimasto a lungo nascosto nella coda poiché i casi comuni non toccavano mai la riga di codice errata.

Una chiamata al foglio di calcolo, cinquanta distribuzioni alla base

Le funzioni coprono le famiglie utilizzate in una cartella di lavoro statistica. Troviamo la famiglia normale, NORM.DIST e NORM.S.DIST con le rispettive inverse; la famiglia gamma e chi-quadrato, GAMMA.DIST, CHISQ.DIST, CHISQ.DIST.RT, CHISQ.INV.RT; la famiglia beta, BETA.DIST e BETA.INV; le distribuzioni campionarie T.DIST, T.DIST.2T, F.DIST e F.INV; la coppia discreta BINOM.DIST e POISSON.DIST; e i supporti di inferenza come CONFIDENCE.T e CONFIDENCE.NORM. Dal punto di vista del chiamante ciascuna rappresenta una singola formula. Si impostano gli input nelle celle, si richiede la valutazione alla cartella di lavoro e si legge il risultato.

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;

Il metodo Calculate della cartella di lavoro compila e valuta una formula ad-hoc rispetto al foglio attivo e restituisce un tipo Variant. Un dettaglio può trarre in inganno al primo tentativo: il parser di formule alla base di Calculate utilizza il punto e virgola come separatore degli argomenti, perciò si scrive =SUM(A1;B1) anziché =SUM(A1,B1). Le formule memorizzate nelle celle conservano la virgola standard di Excel. Lo stesso valutatore gestisce ogni funzione statistica indicata di seguito, perciò il funzionamento di una di esse in Calculate garantisce quello delle altre.

Le due funzioni alla base di tutte le altre

La maggior parte delle distribuzioni cumulative di questo insieme non viene calcolata sommando o integrando le rispettive definizioni. Sono calcolate a partire da due funzioni speciali: la funzione gamma incompleta inferiore regolarizzata, indicata con P(a, x), e la funzione beta incompleta regolarizzata, indicata con I_x(a, b). Internamente questi sono gli helper su cui si appoggiano i dispatcher, e la catena è breve. La CDF del chi-quadrato è la CDF gamma con forma df/2 e scala 2. La CDF gamma corrisponde direttamente a P(a, x). Le funzioni cumulative t, F e binomiale sono tutte valori della beta incompleta regolarizzata ai corretti argomenti. La CDF di Poisson corrisponde alla gamma incompleta superiore Q. Implementando con cura le funzioni gamma e beta, una dozzina di distribuzioni ne eredita l'accuratezza gratuitamente.

Il termine "regolarizzato" costituisce il punto centrale. La funzione gamma incompleta grezza cresce come un fattoriale e l'integrale beta grezzo può subire underflow o overflow molto prima che avvenga per il risultato. Le forme regolarizzate sono divise per la gamma o beta completa, perciò risiedono interamente nell'intervallo da zero a uno, che corrisponde esattamente all'ambito di una probabilità. Tale normalizzazione è ciò che consente alla stessa routine di gestire un chi-quadrato con due gradi di libertà e uno con duecento senza che i termini intermedi superino i limiti di un double. Spiega inoltre perché non si calcola una CDF sommando una lunga coda di termini di densità: ciascun termine porta con sé il proprio errore di arrotondamento, gli errori si accumulano con il procedere della serie, e la funzione speciale regolarizzata evita del tutto la somma valutando invece una serie rapidamente convergente o una frazione continua.

Serie al di sotto della diagonale, frazione continua al di sopra

La routine della gamma incompleta prende una decisione prima di effettuare qualsiasi calcolo: confronta x con a + 1. Tale limite non è arbitrario. Lo sviluppo in serie di potenze di P(a, x) converge rapidamente quando x è piccolo rispetto ad a, e lentamente (fino a diventare inutile) quando x è grande. La frazione continua presenta il comportamento opposto. Pertanto, il motore utilizza lo sviluppo in serie di potenze per x inferiore ad a + 1 e una frazione continua di Lentz per x pari o superiore ad a + 1, assegnando a ciascun ramo solo l'operazione per cui è ottimizzato.

La frazione continua richiede un controllo. Il metodo di Lentz opera gestendo un numeratore e un denominatore parziali e invertendo il denominatore a ogni passaggio; se uno dei due si avvicina a zero, l'inversione fallisce. La soluzione consiste in un limite minimo piccolissimo: ogni volta che un termine intermedio scende al di sotto di circa 1e-30 in modulo, viene bloccato a 1e-30, il che mantiene finita la ricorsione senza alterare il valore convergente. Lo stesso blocco appare nella frazione continua della beta incompleta per il medesimo motivo. Si tratta di una piccola costante che svolge un lavoro fondamentale, facendo la differenza tra una valutazione stabile e una divisione per un valore indistinguibile da zero.

La coda superiore, Q(a, x), corrisponde semplicemente a 1 meno P(a, x), e questo è il modo in cui viene calcolato il ramo cumulativo di Poisson: la probabilità di al massimo k eventi con media λ è pari a Q(k + 1, λ). Ricorrere alla gamma incompleta superiore anziché sommare k + 1 termini di Poisson rappresenta, ancora una volta, la scelta di valutare un'unica espressione convergente invece di accumulare molti piccoli valori.

Masse discrete senza overflow del fattoriale

Le distribuzioni discrete presentano un pericolo diverso. Una massa di probabilità binomiale comporta un coefficiente binomiale, e il coefficiente per 52 su 26 è un intero enorme. Calcolandolo direttamente, il numeratore supera i limiti di un double prima della divisione che lo ricondurrebbe a una probabilità sensata. Il motore non lo calcola mai in questo modo. Calcola i fattoriali nello spazio logaritmico tramite la funzione log-gamma, somma e sottrae i logaritmi, integra il logaritmo delle probabilità di successo e fallimento, ed esegue l'esponenziale una sola volta alla fine.

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

La funzione log-gamma stessa è un'approssimazione di Lanczos, precisa sull'intero asse positivo e rapida da valutare. Poiché ogni valore elevato viene gestito come logaritmo fino all'espressione Exp finale, il valore massimo prodotto dalla routine è la probabilità stessa, che è al massimo pari a uno. La funzione di massa di Poisson segue la stessa logica, con il singolo termine log-gamma a rappresentare il fattoriale al denominatore. Le forme chiuse sono gestite come casi speciali ai limiti, dove p è esattamente zero o uno, in modo che il codice non chiami mai Ln(0). HotXLS restituisce 0.2460938 per BINOM.DIST(5,10,0.5,FALSE) e 0.6766764 per la cumulativa POISSON.DIST(2,2,TRUE), in linea con le cifre mostrate da Excel.

Inverse tramite delimitazione (bracketing) della curva diretta

Una funzione di distribuzione inversa pone la domanda opposta: data una probabilità, trovare il valore la cui CDF corrisponde a essa. Una sola inversa in questo insieme dispone di una formula diretta rapida. NORM.S.INV, l'inversa normale standard, utilizza un'approssimazione razionale di Acklam, una coppia di rapporti polinomiali precisi all'incirca quanto un double sull'intero intervallo, suddivisa in una regione centrale e due code. Si tratta di una valutazione in forma chiusa che non richiede iterazioni.

Le altre inverse non dispongono di tale formula, perciò il motore le inverte numericamente. Delimita la risposta con un limite inferiore e uno superiore scelti dal supporto della distribuzione, quindi dimezza l'intervallo (bisezione): valuta la CDF diretta nel punto medio, sposta il limite che mantiene racchiusa la probabilità target e ripete l'operazione fino a restringere l'intervallo. Per le inverse gamma e chi-quadrato la delimitazione inizia da zero e da una stima superiore ampia basata su forma e scala, raddoppiando il limite superiore se la probabilità non è ancora racchiusa. L'inversa t delimita limiti simmetrici che si allargano verso l'esterno; l'inversa F esegue la bisezione su un intervallo non negativo. Il costo corrisponde a poche decine di valutazioni della CDF per chiamata, impercettibile alla velocità del foglio di calcolo, e il vantaggio risiede nel fatto che ogni inversa ha esattamente la stessa precisione della funzione diretta che inverte. Questo spiega perché un percorso completo come CHISQ.DIST(CHISQ.INV(0.7,5),5,TRUE) restituisca 0.7 con estrema precisione.

Il logaritmo in base 10 nascosto nella coda

Questo è un bug che vale la pena raccontare, poiché è di quelli che sopravvivono a lungo. La routine per la normale inversa di Acklam presenta tre rami. Il ramo centrale ampio, utilizzato quando la probabilità si colloca tra circa 0.025 e 0.975, passa l'input attraverso un rapporto polinomiale che non prevede logaritmi. I due rami delle code, per probabilità molto piccole o molto grandi, calcolano prima un logaritmo dell'input, poiché la coda si comporta come la radice quadrata del logaritmo naturale negativo di p.

Una prima versione del ramo della coda calcolava un logaritmo in base 10 laddove era richiesto il logaritmo naturale. I due differiscono di un fattore costante di circa 2.30, perciò i risultati nella coda erano errati di un margine significativo e costante. Eppure la funzione appariva corretta in ogni verifica superficiale, poiché i controlli comuni si concentrano al centro dell'intervallo. NORM.S.INV(0.5) è pari a zero, NORM.S.INV(0.975) corrisponde al valore classico 1.959964, ed entrambi passano per il polinomio centrale che non richiama alcun logaritmo. L'errore emergeva solo quando una probabilità entrava nella coda, ad esempio con NORM.S.INV(0.001), che deve restituire -3.0902323 e invece forniva un valore distorto dal rapporto tra logaritmo naturale e in base 10. Qualsiasi funzione dipendente dalla normale inversa nella sua coda, inclusi i supporti per gli intervalli di confidenza, ereditava la stessa distorsione. La lezione è semplice ma costosa: una funzione con struttura a rami richiede punti di test all'interno di ciascun ramo, poiché un percorso comune corretto maschera facilmente un ramo secondario errato. La correzione è consistita nel sostituire il logaritmo in base 10 con il logaritmo naturale, allineando i valori della coda a quelli di Excel.

Il segno di x determina la coda della distribuzione t

La funzione cumulativa t di Student presenta una sottigliezza facile da invertire. Il suo valore deriva dalla beta incompleta regolarizzata valutata a df / (df + x²), ma tale valore beta corrisponde alla probabilità nella coda oltre il modulo di x, non alla probabilità cumulativa fino a x. La forma simmetrica della distribuzione t fa sì che la conversione dipenda dal lato dello zero in cui cade 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

Per x superiore a zero la probabilità cumulativa è pari a uno meno la metà della coda simmetrica; per x inferiore a zero corrisponde alla metà di tale coda; a zero è esattamente pari a un mezzo. Restituendo direttamente il valore beta si farebbe riferimento al lato errato della distribuzione, con uno scostamento pari all'intero corpo della curva per ogni x diverso da zero. Le varianti a coda destra e a due code si basano sullo stesso ramo, motivo per cui T.DIST.2T(1,1) restituisce 0.5 e T.DIST(1,1,TRUE) 0.75, e l'inversa T.INV esegue la bisezione rispetto a questa CDF corretta completando il ciclo.

Nulla di tutto questo è visibile dalla cella, ed è questo il risultato desiderato. Scrivi una formula e leggi un numero che concorda con Excel. Se stai estendendo il motore con la tua logica, le procedure per registrare una funzione sono trattate nella nostra guida al motore di formule e alle funzioni personalizzate, e il modo in cui le formule fanno riferimento a diversi fogli e intervalli con nome è descritto nell'articolo sui nomi definiti e le formule tra fogli. Tutto questo è fornito all'interno del componente per fogli di calcolo HotXLS per Delphi e C++Builder, insieme alle API di lettura, scrittura, grafici e formattazione trattate in altre sezioni di questo blog.