Función de distribución acumulativa normal en C/C++

Question

Me preguntaba si había funciones estadísticas incorporadas en las bibliotecas matemáticas que son parte de las bibliotecas estándar de C++ como cmath. Si no, ¿pueden recomendar una buena biblioteca de estadísticas que tenga una función de distribución normal acumulativa? Gracias por adelantado.

Más específicamente, estoy buscando usar/crear una función de distribución acumulativa.

Answer 1

me di cuenta de cómo hacerlo usando GSL, a sugerencia de las personas que respondieron antes que yo, pero luego encontró una solución no-biblioteca (espero que esto ayuda a muchas personas por ahí que busca es como si estuviera):

#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 es tan bueno como el estándar: D aquí va: boost maths/statistical.

Answer 3

Aquí hay una implementación C++ independiente de la distribución normal acumulativa en 14 líneas de código.

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

a partir de muestras de 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 no es una función directa. Pero dado que la función de error de Gauss y su función complementaria se relaciona con la función de distribución acumulativa normal (ver here) podemos usar el c-función implementada erfc:

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

lo uso para los cálculos estadísticos y funciona muy bien. No es necesario usar coeficientes.

Answer 6

Las implementaciones de la CDF normal dada aquí son precisión simple aproximaciones que han tenido float sustituye con double y por lo tanto sólo son precisos a 7 u 8 cifras significativas (decimal).
Para una implementación de VB de doble precisión de Hart aproximación, vea la figura 2 de West Better approximations to cumulative normal functions.

Editar: Mi traducción de aplicación de Occidente en 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; 
} 

Tenga en cuenta que he reordenado expresiones en las formas más familiares para las series y las continuas aproximaciones de fracciones. El último número mágico en el código de West es la raíz cuadrada de 2 π, que he diferido al compilador en la primera línea explotando la identidad acos (0) = & frac12; π.
He comprobado tres veces los números mágicos, pero siempre existe la posibilidad de que haya escrito mal algo. Si detecta un error tipográfico, ¡por favor coméntelo!

Los resultados de los datos de prueba, John Cook, utilizó en su respuesta son

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 

que tomar algún consuelo en el hecho de que están de acuerdo con todos los dígitos dados por los resultados de Mathematica.