Кумулятивная функция нормального распределения в C / C ++

Question

Мне было интересно, есть ли статистические функции, встроенные в математические библиотеки, которые являются частью стандартных библиотек C ++, таких как cmath. Если нет, можете ли вы, ребята, порекомендовать хорошую статистическую библиотеку, в которой была бы кумулятивная функция нормального распределения? Заранее спасибо.

Точнее, я хочу использовать / создать интегральную функцию распределения.

Answer 1

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.

Я использую его для статистических расчетов, и он отлично работает. Нет необходимости использовать коэффициенты.

Answer 2

Вот отдельная реализация 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";
}

Answer 3

Я выяснил, как это сделать, используя 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;
}

Answer 4

Повышение так же хорошо, как стандарт: D здесь вы идете: повышение математики / статистические .

Answer 5

Реализации нормального CDF, приведенные здесь, представляют собой аппроксимации одинарной точности , в которых float заменены на double и, следовательно, точны только до 7 или 8 значащих (десятичных) цифр.
Для реализации VB приближения двойной точности Харта см. Рис. 2 из лучшего приближения Запада к кумулятивным нормальным функциям .

Редактировать : Мой перевод реализации Уэста на 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;
}

Обратите внимание, что я переставил выражения в более привычные формы для приближений рядов и непрерывных дробей. Последнее магическое число в коде Уэста - это квадратный корень из 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.

Answer 6

Из образцов 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;
}

Авторские права 1993-2012 NVIDIA Corporation. Все права защищены.

Answer 7

С 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 в знаменателе для каждой формулы.

Надеюсь, это поможет.