Präzise Berechnung von CDF der Standardnormalverteilung unter Verwendung der C-Standardmath-Bibliothek

Question

Die Standard-C-Math-Bibliothek bietet keine Funktion zur Berechnung der CDF der Standardnormalverteilung, normcdf(). Es bietet jedoch eng verwandte Funktionen: die Fehlerfunktion erf() und die komplementäre Fehlerfunktion erfc(). Der schnellste Weg, um die CDF zu berechnen ist häufig über die Fehlerfunktion, die unter Verwendung von vordefinierten Konstante M_SQRT1_2 darzustellen √½:

double normcdf (double a) 
{ 
    return 0.5 + 0.5 * erf (M_SQRT1_2 * a); 
} 

Offensichtlich Diese leiden unter massiver subtraktiven Löschung in der negativen Halbebene und sind nicht geeignet für die Mehrheit der Anwendungen. Da die Stornierung Problem leicht durch die Verwendung von erfc(), die jedoch in der Regel von etwas Leistung vermieden wird niedriger als erf(), die am häufigsten empfohlene Berechnung ist:

double normcdf (double a) 
{ 
    return 0.5 * erfc (-M_SQRT1_2 * a); 
} 

Einige Tests zeigen, dass die maximale ULP Fehler in der negativen Halb entstanden Flugzeug ist immer noch ziemlich groß. Unter Verwendung einer Doppelpräzisionsimplementierung von erfc(), die auf 0,51 ul genau ist, kann man Fehler von bis zu 1705,44 uls in normcdf() beobachten. Das Problem hierbei ist, dass der Berechnungsfehler in der Eingabe zu erfc() durch die exponentielle Skalierung, die erfc() inhärent ist, vergrößert wird (Siehe diese answer für eine Erklärung der Fehlervergrößerung durch Potenzierung).

Der folgende Beitrag zeigt, wie eine (fast) richtig gerundet, Produkte erreichen können, wenn sie mit beliebiger Genauigkeit Konstanten wie √½ Gleitkommaoperanden Multiplikation:

Nicolas Brisebarre und Jean-Michel Muller, „Korrekt gerundet Multiplikation mit beliebigen Präzisionskonstanten ", IEEE Transactions on Computers, Vol. Die Methode, auf die sich das Papier stützt, beruht auf der fusionierten Multiply-Add-Operation, die bei neueren Implementierungen aller gängigen Prozessorarchitekturen verfügbar ist und in C exponiert ist über die Standard-Mathematikfunktion fma(). Dies führt zu der folgenden Version:

double normcdf (double a) 
{ 
    double SQRT_HALF_HI = 0x1.6a09e667f3bcd0p-01; // 7.0710678118654757e-01 
    double SQRT_HALF_LO = -0x1.bdd3413b264560p-55; // -4.8336466567264567e-17 

    return 0.5 * erfc (fma (-SQRT_HALF_HI, a, -SQRT_HALF_LO * a)); 
} 

Tests zeigen, dass dies den maximalen Fehler im Vergleich zur vorherigen Version ungefähr halbiert. Unter Verwendung der gleichen hochpräzisen erfc() Implementierung wie zuvor ist der maximale beobachtete Fehler 842,71 uls. Dies ist immer noch weit entfernt von dem üblichen Ziel, grundlegende mathematische Funktionen mit einem Fehler von höchstens einigen ul zu versehen.

Gibt es eine effiziente Methode, die die genaue Berechnung von normcdf() ermöglicht, und die nur Funktionen verwendet, die in der C-Standardbibliothek verfügbar sind?

Answer 1

Ein Weg um die Genauigkeitsgrenzen der in der Frage beschriebenen Ansätze zu umgehen, ist die begrenzte Verwendung der Doppel-Doppel-Berechnung. Dies beinhaltet die Berechnung von -sqrt (0.5) * a als ein Paar double Variablen h und l in Kopf/Schwanz-Mode. Der höherwertige Anteil h des Produkts wird an erfc() übergeben, während der niedrigwertige Anteil l dann verwendet wird, um das erfc() Ergebnis basierend auf der lokalen Steigung der komplementären Fehlerfunktion bei h zu interpolieren.

Die Ableitung von erfc (x) ist -2 * exp (-x * x)/√π. Man möchte jedoch die ziemlich teure Berechnung von exp (-x * x) vermeiden.Es ist known, dass für x> 0, erfc (x) ~ = 2 * exp (-x * x)/(√π * (x + sqrt (x * x + 4/π)). Daher asymptotisch erfc '(x) ~ = -2 * x * erfc (x), und daraus folgt, dass für | l | «| h |, erfc (h + 1) ~ = erfc (h) - 2 * h * l * erfc (.. h) die Negation der letztere Begriff leicht in die Berechnung von l gezogen werden kann, gelangt man auf das folgende Implementierung für doppelte Genauigkeit (unter Verwendung von IEEE-754 binary64):

double my_normcdf (double a) 
{ 
    double h, l, r; 
    const double SQRT_HALF_HI = 0x1.6a09e667f3bcd0p-01; // 7.0710678118654757e-01 
    const double SQRT_HALF_LO = -0x1.bdd3413b264560p-55; // -4.8336466567264567e-17 

    /* clamp input as normcdf(x) is either 0 or 1 asymptotically */ 
    if (fabs (a) > 38.625) a = (a < 0.0) ? -38.625 : 38.625; 

    h = fma (-SQRT_HALF_HI, a, -SQRT_HALF_LO * a); 
    l = fma (SQRT_HALF_LO, a, fma (SQRT_HALF_HI, a, h)); 
    r = erfc (h); 
    if (h > 0.0) r = fma (2.0 * h * l, r, r); 
    return 0.5 * r; 
} 

die maximale beobachtete Fehler, die unter Verwendung von gleiche erfc() Implementierung wie zuvor verwendet, ist 1,96 ul. Die entsprechende Single-Precision-Implementierung (unter Verwendung von IEEE-754 binary32) ist:

float my_normcdff (float a) 
{ 
    float h, l, r; 
    const float SQRT_HALF_HI = 0x1.6a09e6p-01f; // 7.07106769e-1 
    const float SQRT_HALF_LO = 0x1.9fcef4p-27f; // 1.21016175e-8 

    /* clamp input as normcdf(x) is either 0 or 1 asymptotically */ 
    if (fabsf (a) > 14.171875f) a = (a < 0.0f) ? -14.171875f : 14.171875f; 

    h = fmaf (-SQRT_HALF_HI, a, -SQRT_HALF_LO * a); 
    l = fmaf (SQRT_HALF_LO, a, fmaf (SQRT_HALF_HI, a, h)); 
    r = erfcf (h); 
    if (h > 0.0f) r = fmaf (2.0f * h * l, r, r); 
    return 0.5f * r; 
}