Excel-Statistikfunktionen in Delphi: NORM, CHISQ, BETA

Geben Sie =NORM.DIST(115,100,15,TRUE) in eine Zelle ein, und Excel gibt ohne Umschweife 0.8413447 zurück. Der Aufruf liest sich wie ein einfacher Tabellen-Lookup. Das ist er aber nicht. Hinter dieser einen Zahl steht die kumulierte Normalverteilung, ein Integral ohne geschlossene Form, und hinter CHISQ.INV.RT und BETA.DIST verbergen sich mathematische Spezialfunktionen, die eine präzise arbeitende Bibliothek berechnen muss, anstatt sie von Hand abzurunden. Eine Tabellenkalkulationskomponente, die Excel-Kompatibilität beansprucht, muss diese Werte bis zur letzten von Excel angezeigten Ziffer reproduzieren, was bedeutet, dass die numerischen Methoden und nicht nur die Funktionsnamen nachgebildet werden müssen.

HotXLS implementiert mehr als fünfzig dieser statistischen Funktionen, und die Arbeit, die sie korrekt macht, ist in der Formelleiste fast vollständig unsichtbar. Dies ist ein Rundgang durch die Funktionsweise der Berechnung durch die Engine: der gemeinsame Spezialfunktionskern, die Verzweigungsentscheidungen, die die Arithmetik stabil halten, und ein Fehler bei der inversen Normalverteilung, der sich lange Zeit im äußeren Bereich (Tail) verbarg, weil der Regelfall die fehlerhafte Zeile niemals ausführte.

Ein Arbeitsblattaufruf, fünfzig Verteilungen dahinter

Die Funktionen decken die Familien ab, die in einer Statistik-Arbeitsmappe benötigt werden. Es gibt die Normalverteilungs-Familie, NORM.DIST und NORM.S.DIST mit ihren Inversen; die Gamma- und Chi-Quadrat-Familie, GAMMA.DIST, CHISQ.DIST, CHISQ.DIST.RT, CHISQ.INV.RT; die Beta-Familie, BETA.DIST und BETA.INV; die Stichprobenverteilungen T.DIST, T.DIST.2T, F.DIST und F.INV; das diskrete Paar BINOM.DIST und POISSON.DIST; sowie die Inferenz-Hilfen wie CONFIDENCE.T und CONFIDENCE.NORM. Aus Sicht des Aufrufers ist jede davon eine einzelne Formel. Sie tragen Werte in Zellen ein, lassen die Arbeitsmappe rechnen und lesen das Ergebnis ab.

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;

Die Methode Calculate der Arbeitsmappe kompiliert und wertet eine Ad-hoc-Formel für das aktive Blatt aus und gibt einen Variant-Wert zurück. Ein Detail bringt Einsteiger beim ersten Versuch oft ins Stolpern: Der Formel-Parser hinter Calculate verwendet das Semikolon als Argumenttrennzeichen, es heißt also =SUM(A1;B1) und nicht =SUM(A1,B1). In Zellen gespeicherte Formeln behalten das von Excel gewohnte Komma bei. Derselbe Evaluator verarbeitet jede unten aufgeführte statistische Funktion. Sobald eine davon in Calculate funktioniert, folgen die restlichen demselben Pfad.

Die zwei Funktionen, auf denen alles andere aufbaut

Die meisten kumulierten Verteilungsfunktionen in diesem Satz werden nicht durch Summieren oder Integrieren ihrer eigenen Definitionen berechnet. Sie werden aus zwei Spezialfunktionen berechnet: der regularisierten unteren unvollständigen Gammafunktion, geschrieben P(a, x), und der regularisierten unvollständigen Betafunktion, geschrieben I_x(a, b). Intern sind dies die Helfer, auf die sich die Dispatcher stützen, und die Verbindungskette ist kurz. Die Chi-Quadrat-Verteilungsfunktion (CDF) ist die Gamma-CDF mit der Form df/2 und der Skalierung 2. Die Gamma-CDF ist direkt P(a, x). Die kumulierten Funktionen von t-, F- und Binomialverteilung sind allesamt Werte der regularisierten unvollständigen Betafunktion bei den passenden Argumenten. Die Poisson-CDF ist die obere unvollständige Gammafunktion Q. Implementiert man die Gamma- und Betafunktion präzise, erben ein Dutzend Verteilungen deren Genauigkeit völlig kostenlos.

Das Wort "regularisiert" trifft den eigentlichen Kern. Die rohe unvollständige Gammafunktion wächst wie eine Fakultät, und das rohe Beta-Integral kann über- oder unterlaufen, lange bevor das Endergebnis dies tut. Die regularisierten Formen sind durch die vollständige Gamma- bzw. Betafunktion dividiert, sodass sie vollständig im Intervall von Null bis Eins liegen – was genau dem Bereich entspricht, den eine Wahrscheinlichkeit einnimmt. Diese Normalisierung sorgt dafür, dass dieselbe Routine eine Chi-Quadrat-Verteilung mit zwei Freiheitsgraden und eine mit zweihundert bedienen kann, ohne dass die Zwischenterme die Grenzen eines Double-Werts überschreiten. Es erklärt auch, warum man eine CDF nicht durch Aufsummieren vieler einzelner Dichtewerte berechnet: Jeder Term bringt seinen eigenen Rundungsfehler mit, die Fehler summieren sich im Verlauf der Reihe auf, und die regularisierte Spezialfunktion umgeht die Summe vollständig, indem sie stattdessen eine schnell konvergierende Reihe oder einen Kettenbruch auswertet.

Reihen unterhalb der Diagonalen, Kettenbruch oberhalb

Die Routine für die unvollständige Gammafunktion trifft eine Entscheidung, bevor sie irgendetwas berechnet: Sie vergleicht x mit a + 1. Diese Grenze ist nicht willkürlich gewählt. Die Potenzreihenentwicklung von P(a, x) konvergiert schnell, wenn x im Verhältnis zu a klein ist, und langsam (schließlich unbrauchbar langsam), wenn x groß ist. Der Kettenbruch verhält sich genau umgekehrt. Die Engine verwendet daher die Potenzreihe für x-Werte unter a + 1 und einen Kettenbruch nach Lentz für x-Werte bei oder über a + 1. So wird jeder Berechnungszweig nur mit der Arbeit betraut, für die er sich am besten eignet.

Der Kettenbruch benötigt eine Schutzmaßnahme. Die Methode nach Lentz arbeitet mit einem laufenden Zähler und Nenner und invertiert den Nenner bei jedem Schritt. Nähert sich einer der beiden Null an, scheitert die Invertierung. Die Lösung ist ein winziger Mindestwert: Wann immer ein Zwischenwert betragsmäßig unter etwa 1e-30 fällt, wird er auf 1e-30 begrenzt, was die Rekursion endlich hält, ohne den konvergierten Wert zu verfälschen. Derselbe Mindestwert wird aus demselben Grund auch im Kettenbruch der unvollständigen Betafunktion angewendet. Es ist eine kleine Konstante mit tragender Rolle – der Unterschied zwischen einer stabilen Auswertung und einer Division durch einen Wert, der sich praktisch nicht von Null unterscheidet.

Der obere Bereich (Tail), Q(a, x), is schlicht 1 minus P(a, x), und so wird auch der kumulierte Zweig für Poisson berechnet: Die Wahrscheinlichkeit von höchsten k Ereignissen bei einem Erwartungswert λ ist Q(k + 1, λ). Der Weg über die obere unvollständige Gammafunktion anstelle des Summierens von k + 1 Poisson-Termen ist wiederum die Entscheidung, einen einzelnen konvergierenden Ausdruck auszuwerten, anstatt viele kleine Werte aufzusummieren.

Diskrete Massen ohne Fakultätsüberlauf

Die diskreten Verteilungen bergen eine andere Gefahr. Eine binomiale Wahrscheinlichkeitsfunktion enthält einen Binomialkoeffizienten, und der Koeffizient für "52 über 26" ist eine enorme Ganzzahl. Bildet man ihn direkt ab, läuft der Zähler bei einem Double-Wert über, noch bevor die Division stattfinden kann, die ihn wieder auf eine sinnvolle Wahrscheinlichkeit zurückführen würde. Die Engine berechnet ihn daher niemals direkt. Sie berechnet die Fakultäten im logarithmischen Raum über die Log-Gamma-Funktion, addiert und subtrahiert die Logarithmen, rechnet die Logarithmen der Erfolgs- und Misserfolgswahrscheinlichkeiten ein und wendet erst ganz am Ende die Exponentialfunktion an.

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

Die Log-Gamma-Funktion selbst ist eine Lanczos-Approximation, die über die gesamte positive Achse hinweg präzise und kostengünstig auszuwerten ist. Da jeder große Wert bis zum abschließenden Exp als Logarithmus gehalten wird, ist die größte Zahl, die die Routine jemals tatsächlich bildet, die Wahrscheinlichkeit selbst (also höchstens eins). Die Poisson-Massenfunktion folgt demselben Rezept, wobei der einzelne Log-Gamma-Term anstelle der Fakultät im Nenner steht. Die geschlossenen Formen werden an den Rändern, wo p genau Null oder Eins ist, als Sonderfälle behandelt, sodass der Code niemals Ln(0) aufruft. HotXLS gibt 0.2460938 für BINOM.DIST(5,10,0.5,FALSE) und 0.6766764 für das kumulierte POISSON.DIST(2,2,TRUE) zurück, was exakt mit den von Excel ausgegebenen Stellen übereinstimmt.

Inverse Verteilungsfunktionen durch Einklammern der Vorwärtskurve

Eine inverse Verteilungsfunktion stellt die umgekehrte Frage: Bestimme für eine gegebene Wahrscheinlichkeit den Wert, dessen kumulierte Verteilungsfunktion (CDF) diesem Wert entspricht. Nur eine Inverse in diesem Satz besitzt eine schnelle direkte Formel. NORM.S.INV, die inverse Standardnormalverteilung, verwendet eine rationale Approximation nach Acklam – ein Paar von Polynomverhältnissen, das über den gesamten Bereich hinweg in etwa der Genauigkeit eines Double-Werts entspricht, aufgeteilt in einen zentralen Bereich und zwei Endbereiche (Tails). Dies ist eine geschlossene Auswertung ohne Iteration.

Die anderen inversen Funktionen haben keine solche Formel, weshalb die Engine sie numerisch invertiert. Sie klammert die Antwort mit einer unteren und oberen Grenze ein, die aus dem Definitionsbereich der Verteilung gewählt werden, und führt dann eine Bisektion (Intervallhalbierung) durch: Wertung der Vorwärts-CDF am Mittelpunkt, Verschieben der Grenze, die die Zielwahrscheinlichkeit eingeschlossen hält, und Wiederholen, bis das Intervall schmal ist. Für die Gamma- und Chi-Quadrat-Inversen beginnt die Einklammerung bei Null und einer großzügigen oberen Schätzung, die auf Form und Skalierung basiert, wobei die obere Grenze verdoppelt wird, wenn die Wahrscheinlichkeit noch nicht eingeschlossen ist. Die inverse t-Verteilung klammert symmetrische Grenzen ein, die sich nach außen erweitern; die inverse F-Verteilung halbiert auf einem nicht-negativen Intervall. Die Kosten belaufen sich auf ein paar Dutzend CDF-Auswertungen pro Aufruf, was bei Tabellenkalkulations-Geschwindigkeit nicht ins Gewicht fällt. Der Vorteil ist jedoch, dass jede Inverse exakt so genau ist wie die Vorwärtsfunktion, die sie invertiert. Aus diesem Grund gibt ein Rundlauf wie CHISQ.DIST(CHISQ.INV(0.7,5),5,TRUE) haargenau 0.7 zurück.

Der Logarithmus zur Basis 10, der sich im Tail verbarg

Eine frühere Version des Zweigs für die Endbereiche verwendete einen Logarithmus zur Basis 10, wo eigentlich der natürliche Logarithmus hingehörte. Die beiden unterscheiden sich um einen konstanten Faktor von etwa 2,30, weshalb die Ergebnisse im Tail-Bereich um eine gleichbleibende, erhebliche Abweichung falsch waren. Und dennoch sah die Funktion bei jeder oberflächlichen Prüfung in Ordnung aus, weil solche Prüfungen im mittleren Bereich stattfinden. NORM.S.INV(0.5) ist Null, NORM.S.INV(0.975) entspricht den standardmäßigen 1.959964, und beide durchlaufen das zentrale Polynom, das überhaupt keinen Logarithmus aufruft. Der Fehler trat erst auf, wenn eine Wahrscheinlichkeit in den Endbereich rutschte – etwa bei NORM.S.INV(0.001), was eigentlich -3.0902323 zurückgeben müsste, aber stattdessen einen durch das Verhältnis von natürlichem Logarithmus zu Logarithmus zur Basis 10 verfälschten Wert lieferte. Jede Funktion, die in ihren Endbereichen von der inversen Normalverteilung abhängt (einschließlich der Konfidenzintervall-Helfer), erbte dieselbe Abweichung. Die Lehre ist banal und kostspielig: Eine Funktion mit Verzweigungsstruktur benötigt Testpunkte in jedem einzelnen Zweig, da ein korrekter Hauptpfad einen fehlerhaften Randpfad problemlos kaschieren kann. Die Behebung war eine Änderung von einem einzigen Token – vom Logarithmus zur Basis 10 zum natürlichen Logarithmus –, woraufhin die Werte im Endbereich exakt auf die Werte von Excel einschnappten.

Das Vorzeichen von x entscheidet über den Tail-Bereich der t-Verteilung

Die kumulierte Student-t-Verteilungsfunktion enthält eine Feinheit, die man leicht verwechseln kann. Ihr Wert stammt aus der regularisierten unvollständigen Betafunktion, ausgewertet bei df / (df + x²). Dieser Beta-Wert entspricht jedoch der Wahrscheinlichkeit im Endbereich jenseits der Größe von x, nicht der kumulierten Wahrscheinlichkeit bis zu x. Die symmetrische Form der t-Verteilung bedeutet, dass die Umwandlung davon abhängt, auf welcher Seite von Null x liegt.

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

Für x-Werte über Null ist die kumulierte Wahrscheinlichkeit eins minus der halbe symmetrische Endbereich; für x-Werte unter Null ist sie die Hälfte dieses Endbereichs; bei Null ist sie genau ein Halb. Gibt man den Beta-Wert direkt zurück, meldet man die falsche Seite der Verteilung – für jedes x ungleich Null um den gesamten Körper der Kurve verschoben. Die rechtsseitigen und zweiseitigen Varianten bauen auf demselben Zweig auf, weshalb T.DIST.2T(1,1) als 0.5 und T.DIST(1,1,TRUE) als 0.75 zurückgegeben wird, und die Inverse T.INV führt die Bisektion gegen diese korrigierte CDF durch, sodass der Rundlauf geschlossen wird.

Nichts davon ist in der Zelle sichtbar, und das ist das gewünschte Ergebnis. Sie schreiben eine Formel und lesen eine Zahl ab, die mit Excel übereinstimmt. Wenn Sie die Engine mit Ihrer eigenen Logik erweitern möchten, werden die Mechanismen zur Registrierung einer Funktion in unserer Einführung in die Formel-Engine und benutzerdefinierte Funktionen behandelt, und der Weg, wie Formeln über Tabellen und benannte Bereiche hinweg aufgelöst werden, wird im Artikel über definierte Namen und tabellenübergreifende Formeln beschrieben. All dies wird in der HotXLS-Tabellenkalkulationskomponente für Delphi und C++Builder ausgeliefert, zusammen mit den Lese-, Schreib-, Diagramm- und Formatierungs-APIs, die an anderer Stelle in diesem Blog behandelt werden.