c++ c math statistics distribution

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



math statistics (6)

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.

A partir de muestras 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. Todos los derechos reservados.

Aquí hay una implementación independiente de C ++ 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"; }

Descubrí cómo hacerlo usando gsl, por sugerencia de las personas que respondieron antes que yo, pero luego encontré una solución que no es de biblioteca (con suerte, esto ayuda a muchas personas que lo buscan como yo):

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

Las implementaciones de la CDF normal dadas aquí son aproximaciones de precisión simple que han tenido float reemplazada por double y por lo tanto solo son precisas para 7 u 8 cifras significativas (decimales).
Para una implementación de VB de la aproximación de doble precisión de Hart, vea la figura 2 de las mejores aproximaciones de West a las funciones normales acumulativas .

Editar : Mi traducción de la implementación de West 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 reorganizado expresiones en las formas más familiares para series y aproximaciones de fracciones continuas. 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) = ½ π.
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 que usó John Cook 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

Me reconforta un poco el hecho de que están de acuerdo con todos los dígitos dados para los resultados de Mathematica.

No hay una función directa. Pero dado que la función de error gaussiano y su función complementaria está relacionada con la función de distribución acumulativa normal (ver here ), podemos usar la función c ejecutada erfc :

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

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