Fonction de Distribution normale cumulée en C/C++
je me demandais s'il y avait des fonctions statistiques intégrées dans les bibliothèques mathématiques qui font partie des bibliothèques C++ standard comme cmath. Si ce n'est pas le cas, Pouvez-vous recommander une bonne bibliothèque de statistiques qui aurait une fonction de distribution normale cumulative? Merci à l'avance.
Plus précisément, je cherche à utiliser/créer une fonction de distribution cumulée.
6 réponses
Voici une implémentation autonome C++ de la distribution normale cumulative en 14 lignes de code.
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";
}
Theres n'est pas une fonction droite. Mais puisque la fonction d'erreur gaussienne et sa fonction complémentaire est liée à la fonction de distribution cumulative normale (voir ici) nous pouvons utiliser la fonction C implémentée erfc:
double normalCFD(double value)
{
return 0.5 * erfc(-value * M_SQRT1_2);
}
Je l'utilise pour des calculs statistiques et ça marche très bien. Pas besoin d'utiliser des coefficients.
Boost est aussi bon que le standard :D ici, vous allez: booster de mathématiques/statistiques.
j'ai trouvé comment le faire en utilisant gsl, à la suggestion des gens qui ont répondu avant moi, mais a ensuite trouvé une solution non-bibliothèque (espérons que cela aide beaucoup de gens là-bas qui sont à la recherche pour lui comme je l'étais):
#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;
}
les implémentations du CDF normal donné ici sont précision approximations qui ont eu float remplacé par double et ne sont donc exacts qu'à 7 ou 8 chiffres significatifs (décimaux).
Pour une implémentation VB de Hart's double précision rapprochement, voir la figure 2 de l'Ouest meilleures approximations des fonctions normales cumulatives.
Modifier: ma traduction de West's l'implémentation 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;
}
notez que j'ai réarrangé des expressions dans des formes plus familières pour des approximations de séries et de fractions continues. Le dernier nombre magique dans le code de West est la racine carrée de 2π, que j'ai reportée au compilateur sur la première ligne en exploitant l'identité acos(0) = ½ π.
J'ai vérifié trois fois les nombres magiques, mais il y a toujours une chance que je me sois trompé. Si vous repérez une faute de frappe, merci commentaire!
Les résultats pour les données de test John Cook utilisé dans sa réponse sont
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
je prends un peu de confort du fait qu'ils acceptent tous les chiffres donnés pour les résultats de Mathematica.
From NVIDIA CUDA samples:
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. Tous droits réservés.