Мне было интересно, есть ли статистические функции, встроенные в математические библиотеки, которые являются частью стандартных библиотек С++, таких как cmath. Если нет, можете ли вы, ребята, рекомендовать хорошую статистику, которая будет иметь кумулятивную нормальную функцию распределения? Спасибо заранее.
В частности, я хочу использовать/создать кумулятивную функцию распределения.
Я понял, как это сделать, используя gsl, по предложению людей, которые отвечали передо мной, но затем нашел небиблиотечное решение (надеюсь, это помогает многим людям, которые ищут его, как я):
#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;
}
Theres не прямая функция. Но поскольку гауссовская функция ошибок и ее дополнительная функция связаны с нормальной кумулятивной функцией распределения (см. Здесь или здесь), мы можем использовать реализованную c-функцию erfc (дополнительная функция ошибок):
double normalCDF(double value)
{
return 0.5 * erfc(-value * M_SQRT1_2);
}
Который рассматривает отношение erfc(x) = 1-erf(x) с M_SQRT1_2= √0,5.
Я использую его для статистических расчетов, и он прекрасно работает. Нет необходимости использовать коэффициенты.
Здесь автономная C++ реализация кумулятивного нормального распределения в 14 строках кода.
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";
}
Boost так же хорош, как и стандарт: D здесь вы идете: повысить математику/статистику.
Реализации нормального CDF, приведенного здесь, являются одноточечными приближениями, которые <<20 > заменены на double и, следовательно, являются точными только для 7 или 8 значимых (десятичных) цифр.
Для VB-реализации приближения с двойной точностью Hart см. Рис. 2 "Запад" Лучшее приближение к кумулятивным нормальным функциям.
Изменить. Мой перевод реализации Запада на С++:
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;
}
Обратите внимание, что у меня есть переупорядоченные выражения в более знакомые формы для рядов и продолжения дробных приближений. Последнее магическое число в западном коде - это квадратный корень из 2 & pi;, который я отложил к компилятору в первой строке, используя идентификаторы acos (0) = & frac12; & Пи;.
Я тройной проверил магические числа, но всегда есть вероятность, что я что-то опечалил. Если вы заметили опечатку, прокомментируйте!
Результаты тестовых данных, которые Джон Кук использовал в своем ответе,
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
Я немного успокаиваюсь из-за того, что они согласны со всеми цифрами, данными для результатов Mathematica.
Из образцов 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;
}
С https://en.cppreference.com/w/cpp/numeric/math/erfc
Нормальный CDF можно рассчитать, как показано ниже:
#include <iostream>
#include <cmath>
#include <iomanip>
using namespace std;
double normalCDF(double x) // Phi(-∞, x) aka N(x)
{
return erfc(-x / sqrt(2))/2;
}
Дополнительное примечание, из моего предыдущего опыта, если у вас есть проблема с получением только целых чисел вместо десятичных дробей, использование 2.0 вместо 2 в знаменателе для каждой формулы поможет.
Надеюсь, это поможет.