Fonctions statistiques Excel dans Delphi : NORM, CHISQ, BETA

Saisissez =NORM.DIST(115,100,15,TRUE) dans une cellule et Excel renvoie 0.8413447 sans autre formalité. L'appel ressemble à une simple recherche dans une table. Ce n'est pas le cas. Derrière ce nombre unique se cache la distribution normale cumulative, une intégrale sans forme close, et derrière CHISQ.INV.RT et BETA.DIST résident des fonctions spéciales qu'une bibliothèque rigoureuse doit évaluer, et non estimer à la main. Un composant de feuille de calcul qui revendique la compatibilité avec Excel doit reproduire ces valeurs jusqu'au dernier chiffre affiché par Excel, ce qui implique de reproduire les méthodes numériques, et pas seulement les noms de fonctions.

HotXLS implémente plus de cinquante de ces fonctions statistiques, et le travail qui garantit leur exactitude est presque entièrement invisible depuis la barre de formule. Voici un aperçu de la manière dont le moteur les calcule : le cœur de fonctions spéciales partagé, les décisions de branchement qui maintiennent la stabilité de l'arithmétique, et un bogue de loi normale inverse qui est resté longtemps caché dans la queue (tail) car le cas général ne sollicitait jamais la ligne défectueuse.

Un seul appel de feuille de calcul, cinquante distributions derrière

Les fonctions couvrent les familles qu'un classeur de statistiques sollicite. On y trouve la famille normale (NORM.DIST, NORM.S.DIST et leurs inverses), la famille gamma et chi-carré (GAMMA.DIST, CHISQ.DIST, CHISQ.DIST.RT, CHISQ.INV.RT), la famille bêta (BETA.DIST et BETA.INV), les distributions d'échantillonnage (T.DIST, T.DIST.2T, F.DIST et F.INV), la paire discrète (BINOM.DIST et POISSON.DIST) et les assistants d'inférence tels que CONFIDENCE.T et CONFIDENCE.NORM. Du point de vue de l'appelant, chacune est une formule unique. Vous définissez les entrées dans les cellules, demandez au classeur d'évaluer et lisez le résultat.

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;

La méthode Calculate du classeur compile et évalue une formule ad hoc par rapport à la feuille active et renvoie un type Variant. Un détail surprend lors du premier essai : l'analyseur de formule derrière Calculate prend le point-virgule comme séparateur d'arguments, il s'écrit donc sous la forme =SUM(A1;B1), et non =SUM(A1,B1). Les formules de cellules stockées conservent la virgule standard d'Excel. Le même évaluateur distribue chaque fonction statistique ci-dessous, ainsi, dès que l'une d'elles fonctionne dans Calculate, les autres suivent le même chemin.

Les deux fonctions sur lesquelles repose tout le reste

La plupart des distributions cumulatives de cet ensemble ne sont pas calculées en sommant ou en intégrant leurs propres définitions. Elles sont calculées à partir de deux fonctions spéciales : la fonction gamma incomplète inférieure régularisée, notée P(a, x), et la fonction bêta incomplète régularisée, notée I_x(a, b). En interne, ce sont les assistants sur lesquels s'appuient les répartiteurs, et la chaîne est courte. La fonction de distribution cumulative (CDF) du chi-carré est la CDF gamma avec le paramètre de forme df/2 et d'échelle 2. La CDF gamma est directement P(a, x). Les fonctions cumulatives de Student (t), de Fisher (F) et binomiale sont toutes des valeurs de la fonction bêta incomplète régularisée évaluée avec les bons arguments. La CDF de Poisson est la fonction gamma incomplète supérieure Q. Implémentez correctement les fonctions gamma et bêta, et une douzaine de distributions hériteront gratuitement de leur précision.

Le terme « régularisé » résume tout l'intérêt. La fonction gamma incomplète brute croît comme une factorielle et l'intégrale bêta brute peut s'annuler par sous-passement (underflow) ou dépasser sa capacité (overflow) bien avant le résultat final. Les formes régularisées sont divisées par la fonction gamma ou bêta complète, de sorte qu'elles se situent entièrement dans l'intervalle de zéro à un, ce qui correspond exactement à la plage d'une probabilité. Cette normalisation permet à la même routine de traiter un chi-carré à deux degrés de liberté et un autre à deux cents sans que les termes intermédiaires ne dépassent la limite d'un type double. Elle explique également pourquoi on ne calcule pas une CDF en additionnant une longue queue de termes de densité : chaque terme porte sa propre erreur d'arrondi, les erreurs s'accumulent au fil de la série, et la fonction spéciale régularisée évite totalement cette somme en évaluant à la place une série rapidement convergente ou une fraction continue.

Série sous la diagonale, fraction continue au-dessus

La routine de la fonction gamma incomplète prend une décision avant tout calcul : elle compare x à a + 1. Cette frontière n'est pas arbitraire. Le développement en série entière de P(a, x) converge rapidement lorsque x est petit par rapport à a, et lentement, puis inutilement, lorsque x est grand. La fraction continue présente le caractère inverse. Le moteur utilise donc le développement en série pour x inférieur à a + 1 et une fraction continue de Lentz pour x supérieur ou égal à a + 1, chaque branche n'effectuant que le travail pour lequel elle est efficace.

La fraction continue nécessite une protection. La méthode de Lentz fonctionne en conservant un numérateur et un dénominateur courants et en inversant le dénominateur à chaque étape ; si l'un d'eux s'approche de zéro, l'inversion échoue. La solution consiste en un seuil minimal : chaque fois qu'un terme intermédiaire tombe en dessous d'environ 1e-30 en valeur absolue, il est fixé à 1e-30, ce qui maintient la récurrence finie sans perturber la valeur convergée. La même limitation apparaît dans la fraction continue de la fonction bêta incomplète pour la même raison. C'est une constante minuscule réalisant un travail crucial, faisant la différence entre une évaluation stable et une division par une valeur indiscernable de zéro.

La queue supérieure, Q(a, x), vaut simplement 1 moins P(a, x), et c'est ainsi que la branche cumulative de Poisson est calculée : la probabilité d'avoir au plus k événements de moyenne λ est Q(k + 1, λ). Choisir de passer par la fonction gamma incomplète supérieure plutôt que de sommer k + 1 termes de Poisson revient, là encore, à évaluer une expression convergente au lieu d'accumuler de nombreux petits termes.

Masses discrètes sans dépassement de factorielle

Les distributions discrètes soulèvent un danger différent. Une fonction de masse de probabilité binomiale implique un coefficient binomial, et le coefficient pour 26 parmi 52 est un entier gigantesque. Si on le calculait directement, le numérateur dépasserait la capacité d'un type double avant la division censée le ramener à une probabilité raisonnable. Le moteur ne le calcule jamais directement. Il calcule les factorielles dans l'espace logarithmique via la fonction log-gamma, additionne et soustrait les logarithmes, intègre le logarithme des probabilités de succès et d'échec, et applique l'exponentielle une seule fois à la toute fin.

// 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 fonction log-gamma elle-même est une approximation de Lanczos, précise sur tout l'axe positif et peu coûteuse à évaluer. Comme chaque grande quantité est conservée sous forme de logarithme jusqu'au Exp final, le plus grand nombre que la routine matérialise est la probabilité elle-même, qui vaut au plus un. La fonction de masse de Poisson suit la même recette, avec l'unique terme log-gamma remplaçant la factorielle au dénominateur. Les formes closes font l'objet de cas particuliers aux limites, là où p vaut exactement zéro ou un, afin que le code n'appelle jamais Ln(0). HotXLS renvoie 0.2460938 pour BINOM.DIST(5,10,0.5,FALSE) et 0.6766764 pour la formule cumulative POISSON.DIST(2,2,TRUE), correspondant aux chiffres affichés par Excel.

Inverses par encadrement de la courbe directe

Une fonction de distribution inverse pose la question inverse : étant donné une probabilité, trouver la valeur dont la CDF lui correspond. Seule une fonction inverse dans cet ensemble possède une formule directe rapide. NORM.S.INV, la loi normale standard inverse, utilise une approximation rationnelle d'Acklam, à savoir une paire de rapports polynomiaux précis à peu près à la précision d'un double sur toute la plage, répartie entre une zone centrale et deux queues. Il s'agit d'une évaluation sous forme close sans itération.

Les autres fonctions inverses ne disposant pas d'une telle formule, le moteur procède par inversion numérique. Il encadre la réponse avec une limite inférieure et supérieure choisie dans le support de la distribution, puis procède par dichotomie : il évalue la CDF directe au point médian, déplace la limite qui maintient la probabilité cible encadrée, et répète l'opération jusqu'à ce que l'intervalle soit étroit. Pour les inverses du gamma et du chi-carré, l'encadrement commence à zéro et une estimation supérieure généreuse construite à partir de la forme et de l'échelle, doublant la limite supérieure si la probabilité n'est pas encore encadrée. L'inverse de Student (t) encadre des limites symétriques qui s'élargissent vers l'extérieur ; l'inverse de Fisher (F) procède par dichotomie sur un intervalle non négatif. Le coût correspond à quelques dizaines d'évaluations de CDF par appel, ce qui reste invisible à la vitesse d'exécution d'une feuille de calcul, et l'avantage est que chaque inverse est exactement aussi précis que la fonction directe qu'il inverse. C'est pourquoi un aller-retour tel que CHISQ.DIST(CHISQ.INV(0.7,5),5,TRUE) renvoie 0.7 à un cheveu près.

Le logarithme décimal qui se cachait dans la queue

Une première version de la branche de queue utilisait un logarithme décimal (base 10) là où le logarithme népérien était requis. Les deux diffèrent d'un facteur constant d'environ 2.30, de sorte que les résultats dans les queues étaient erronés d'une marge constante et notable. Pourtant, la fonction semblait correcte lors de chaque vérification informelle, car celles-ci se situent au milieu. NORM.S.INV(0.5) vaut zéro, NORM.S.INV(0.975) vaut la valeur classique 1.959964, et ces deux calculs passent par le polynôme central qui n'appelle aucun logarithme. L'erreur n'apparaissait que lorsqu'une probabilité entrait dans une queue, par exemple NORM.S.INV(0.001), qui doit renvoyer -3.0902323 et renvoyait à la place une valeur décalée par le rapport entre logarithme népérien et décimal. Toute fonction dépendant de la loi normale inverse dans sa queue, y compris les assistants d'intervalle de confiance, héritait du même biais. La leçon est banale et coûteuse : une fonction dotée d'une structure de branchement nécessite des points de test dans chaque branche, car un chemin commun correct masquera volontiers une branche défectueuse. La correction a consisté à remplacer le logarithme en base 10 par le logarithme népérien, et les valeurs de queue se sont alignées sur celles d'Excel.

Le signe de x décide de la queue de la distribution de Student (t)

La fonction cumulative de Student (t) comporte une subtilité facile à inverser. Sa valeur provient de la fonction bêta incomplète régularisée évaluée à df / (df + x²), mais cette valeur bêta représente la probabilité dans la queue au-delà de la valeur absolue de x, et non la probabilité cumulative jusqu'à x. La forme symétrique de la distribution de Student implique que la conversion dépend du côté de zéro où se situe 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

Pour un x supérieur à zéro, la probabilité cumulative vaut un moins la moitié de la queue symétrique ; pour un x inférieur à zéro, elle vaut la moitié de cette queue ; à zéro, elle vaut exactement la moitié. Renvoyez directement la valeur bêta et vous obtiendrez le mauvais côté de la distribution, décalé de tout le corps de la courbe pour tout x non nul. Les variantes à queue droite et à deux queues reposent sur la même branche, c'est pourquoi T.DIST.2T(1,1) renvoie 0.5 et T.DIST(1,1,TRUE) renvoie 0.75, et l'inverse T.INV effectue sa dichotomie par rapport à cette CDF corrigée pour que l'aller-retour concorde.

Rien de tout cela n'est visible depuis la cellule, et c'est le résultat escompté. Vous écrivez une formule et lisez un nombre conforme à Excel. Si vous étendez le moteur avec votre propre logique, les mécanismes d'enregistrement d'une fonction sont traités dans notre guide du moteur de formule et des fonctions personnalisées, et la manière dont les formules ciblent d'autres feuilles et plages nommées est présentée dans l'article sur les noms définis et les formules multi-feuilles. Tout cela est fourni au sein du composant de feuille de calcul HotXLS pour Delphi et C++Builder, avec les API de lecture, d'écriture, de création de graphiques et de formatage traitées par ailleurs sur ce blog.