webentwicklung-frage-antwort-db.com.de

Kumulative Normalverteilungsfunktion in C / C ++

Ich habe mich gefragt, ob Statistikfunktionen in Mathematikbibliotheken integriert sind, die zu den Standard-C++ - Bibliotheken wie cmath gehören. Wenn nicht, könnt ihr eine gute Statistikbibliothek empfehlen, die eine kumulative Normalverteilungsfunktion hat? Danke im Voraus.

Insbesondere möchte ich eine kumulative Verteilungsfunktion verwenden/erstellen.

c++cmathstatisticsdistribution
46
Tyler Brock

Auf Anregung der Leute, die vor mir geantwortet haben, habe ich herausgefunden, wie es mit gsl geht, aber dann habe ich eine Lösung ohne Bibliothek gefunden (hoffentlich hilft das vielen Leuten da draußen, die danach suchen, wie ich es war):

#ifndef Pi 
#define Pi 3.141592653589793238462643 
#endif 

double cnd_manual(double x)
{
  double L, K, w ;
  /* constants */
  double const a1 = 0.31938153, a2 = -0.356563782, a3 = 1.781477937;
  double const a4 = -1.821255978, a5 = 1.330274429;

  L = fabs(x);
  K = 1.0 / (1.0 + 0.2316419 * L);
  w = 1.0 - 1.0 / sqrt(2 * Pi) * exp(-L *L / 2) * (a1 * K + a2 * K *K + a3 * pow(K,3) + a4 * pow(K,4) + a5 * pow(K,5));

  if (x < 0 ){
    w= 1.0 - w;
  }
  return w;
}
9
Tyler Brock

Theres ist keine gerade Funktion. Da jedoch die Gaußsche Fehlerfunktion und ihre Komplementärfunktion mit der normalen kumulativen Verteilungsfunktion verwandt sind (siehe hier oder hier ), können wir die implementierte c-Funktion erfc (komplementäre Fehlerfunktion):

double normalCDF(double value)
{
   return 0.5 * erfc(-value * M_SQRT1_2);
}

Was die Beziehung von erfc(x) = 1-erf(x) mit M_SQRT1_2 = √0,5 berücksichtigt.

Ich benutze es für statistische Berechnungen und es funktioniert großartig. Koeffizienten müssen nicht verwendet werden.

36
JFS

Hier ist eine eigenständige C++ - Implementierung der kumulativen Normalverteilung in 14 Codezeilen.

http://www.johndcook.com/cpp_phi.html

#include <cmath>

double phi(double x)
{
    // constants
    double a1 =  0.254829592;
    double a2 = -0.284496736;
    double a3 =  1.421413741;
    double a4 = -1.453152027;
    double a5 =  1.061405429;
    double p  =  0.3275911;

    // Save the sign of x
    int sign = 1;
    if (x < 0)
        sign = -1;
    x = fabs(x)/sqrt(2.0);

    // A&S formula 7.1.26
    double t = 1.0/(1.0 + p*x);
    double y = 1.0 - (((((a5*t + a4)*t) + a3)*t + a2)*t + a1)*t*exp(-x*x);

    return 0.5*(1.0 + sign*y);
}

void testPhi()
{
    // Select a few input values
    double x[] = 
    {
        -3, 
        -1, 
        0.0, 
        0.5, 
        2.1 
    };

    // Output computed by Mathematica
    // y = Phi[x]
    double y[] = 
    { 
        0.00134989803163, 
        0.158655253931, 
        0.5, 
        0.691462461274, 
        0.982135579437 
    };

        int numTests = sizeof(x)/sizeof(double);

    double maxError = 0.0;
    for (int i = 0; i < numTests; ++i)
    {
        double error = fabs(y[i] - phi(x[i]));
        if (error > maxError)
            maxError = error;
    }

        std::cout << "Maximum error: " << maxError << "\n";
}
32
John D. Cook

Boost ist so gut wie der Standard: D Los geht's: Boost Mathe/Statistik .

9
Hassan Syed

Die hier angegebenen Implementierungen der normalen CDF sind einfache Genauigkeit Näherungen, bei denen float durch double ersetzt wurde und die daher nur auf 7 oder 8 signifikante Stellen (Dezimalstellen) genau sind ) Zahlen.
Für eine VB Implementierung von Harts doppelte Genauigkeit Näherung siehe Abbildung 2 von Wests Bessere Näherungen an kumulative Normalfunktionen =.

Bearbeiten : Meine Übersetzung der Implementierung von West in C++:

double
phi(double x)
{
  static const double RT2PI = sqrt(4.0*acos(0.0));

  static const double SPLIT = 7.07106781186547;

  static const double N0 = 220.206867912376;
  static const double N1 = 221.213596169931;
  static const double N2 = 112.079291497871;
  static const double N3 = 33.912866078383;
  static const double N4 = 6.37396220353165;
  static const double N5 = 0.700383064443688;
  static const double N6 = 3.52624965998911e-02;
  static const double M0 = 440.413735824752;
  static const double M1 = 793.826512519948;
  static const double M2 = 637.333633378831;
  static const double M3 = 296.564248779674;
  static const double M4 = 86.7807322029461;
  static const double M5 = 16.064177579207;
  static const double M6 = 1.75566716318264;
  static const double M7 = 8.83883476483184e-02;

  const double z = fabs(x);
  double c = 0.0;

  if(z<=37.0)
  {
    const double e = exp(-z*z/2.0);
    if(z<SPLIT)
    {
      const double n = (((((N6*z + N5)*z + N4)*z + N3)*z + N2)*z + N1)*z + N0;
      const double d = ((((((M7*z + M6)*z + M5)*z + M4)*z + M3)*z + M2)*z + M1)*z + M0;
      c = e*n/d;
    }
    else
    {
      const double f = z + 1.0/(z + 2.0/(z + 3.0/(z + 4.0/(z + 13.0/20.0))));
      c = e/(RT2PI*f);
    }
  }
  return x<=0.0 ? c : 1-c;
}

Beachten Sie, dass ich Ausdrücke in die bekannteren Formen für Serien und fortlaufende Brüche umgeordnet habe. Die letzte magische Zahl in Wests Code ist die Quadratwurzel von 2π, die ich durch Ausnutzung der Identität acos (0) = ½ π an den Compiler in der ersten Zeile weitergeleitet habe.
Ich habe die magischen Zahlen dreimal überprüft, aber es besteht immer die Möglichkeit, dass ich etwas falsch geschrieben habe. Wenn Sie einen Tippfehler entdecken, kommentieren Sie ihn bitte!

Die Ergebnisse für die Testdaten, die John Cook in seiner Antwort verwendete, sind

 x               phi                Mathematica
-3     1.3498980316301150e-003    0.00134989803163
-1     1.5865525393145702e-001    0.158655253931
 0     5.0000000000000000e-001    0.5
0.5    6.9146246127401301e-001    0.691462461274
2.1    9.8213557943718344e-001    0.982135579437

Ich tröste mich ein wenig damit, dass sie mit allen Ziffern übereinstimmen, die für die Mathematica-Ergebnisse angegeben wurden.

8
thus spake a.k.

Aus NVIDIA CUDA-Beispielen:

static double CND(double d)
{
    const double       A1 = 0.31938153;
    const double       A2 = -0.356563782;
    const double       A3 = 1.781477937;
    const double       A4 = -1.821255978;
    const double       A5 = 1.330274429;
    const double RSQRT2PI = 0.39894228040143267793994605993438;

    double
    K = 1.0 / (1.0 + 0.2316419 * fabs(d));

    double
    cnd = RSQRT2PI * exp(- 0.5 * d * d) *
          (K * (A1 + K * (A2 + K * (A3 + K * (A4 + K * A5)))));

    if (d > 0)
        cnd = 1.0 - cnd;

    return cnd;
}

Copyright 1993-2012 NVIDIA Corporation. Alle Rechte vorbehalten.

5
serbaut

Von https://en.cppreference.com/w/cpp/numeric/math/erfc

Normaler CDF kann wie folgt berechnet werden:

#include <iostream>
#include <cmath>
#include <iomanip>
using namespace std;

double normalCDF(double x) // Phi(-∞, x) aka N(x)
{
    return erfc(-x / sqrt(2))/2;
}

Eine Randnotiz aus meiner bisherigen Erfahrung, wenn Sie ein Problem damit haben, nur Ganzzahlen anstelle von Dezimalzahlen zu erhalten, hilft es, 2,0 anstelle von 2 im Nenner für jede Formel zu verwenden.

Ich hoffe, das hilft.

0
Manohar Reddy Poreddy