Funzione di distribuzione normale cumulativa in C/C++

Question

Mi chiedevo se esistessero funzioni statistiche incorporate in librerie matematiche che fanno parte delle librerie C++ standard come il cmath. In caso contrario, puoi raccomandare una buona libreria di statistiche che avrebbe una funzione di distribuzione normale cumulativa? Grazie in anticipo.

In particolare, sto cercando di utilizzare/creare una funzione di distribuzione cumulativa.

Answer 1

ho capito come farlo usando GSL, su suggerimento delle persone che hanno risposto prima di me, ma poi trovato una soluzione non-library (spero che questo aiuta molte persone là fuori che sono alla ricerca di esso come se fossi):

#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; 
} 

Answer 2

Boost è buono come lo standard: D qui vai: boost maths/statistical.

Answer 3

Ecco un'implementazione C++ stand-alone della distribuzione normale cumulativa in 14 righe di codice.

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"; 
} 

Answer 4

Da campioni di NVIDIA CUDA:

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. All rights reserved.

Answer 5

Theres non ha una funzione diritta. Ma poiché la funzione di errore gaussiana e la sua funzione complementare è legato alla normale funzione di distribuzione cumulativa (vedi here) possiamo usare il c-funzione implementata erfc:

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

lo uso per i calcoli statistici e funziona grande. Non c'è bisogno di usare coefficienti.

Answer 6

Le implementazioni del normale CDF qui proposta sono precisione singola approssimazioni che hanno avuto float sostituito con double e quindi sono solo precisi a 7 o 8 cifre significative (decimale).
Per un'implementazione VB di Hart doppia precisione, vedere la figura 2 di West Better approximations to cumulative normal functions.

Edit: La mia traduzione di implementazione Ovest 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; 
} 

Nota che ho riarrangiato espressioni nelle forme più familiari per le serie e approssimazioni frazione continua. L'ultimo numero magico nel codice di West è la radice quadrata di 2 π, che ho rimandato al compilatore sulla prima riga sfruttando l'identità acos (0) = & frac12; π.
Ho triplicato controllato i numeri magici, ma c'è sempre la possibilità che abbia digitato male qualcosa. Se noti un errore di battitura, per favore commenta!

I risultati per i dati del test John Cook utilizzato nella sua risposta sono

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 

prendo qualche piccolo conforto dal fatto di essere d'accordo a tutte le cifre fornite per i risultati di Mathematica.