C/C++에서의 누적 정규 분포 함수

C/C++에서의 누적 정규 분포 함수

cmath와 같은 표준 C++ 라이브러리의 일부인 수학 라이브러리에 통계 기능이 내장되어 있는지 궁금합니다.그렇지 않다면, 누적 정규 분포 함수를 가진 좋은 통계 라이브러리를 추천해 주시겠습니까?잘 부탁드립니다.

구체적으로는 누적분포함수를 사용/작성하려고 합니다.

이 함수는 직선 함수가 아닙니다.그러나 가우스 오차 함수와 그 보완 함수는 정규 누적 분포 함수와 관련이 있으므로(여기 또는 여기 참조) 구현된 c-함수를 사용할 수 있습니다.erfc(보완적 오류 함수):

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

이 부분은 다음 사항의 관계를 고려합니다.erfc(x) = 1-erf(x)와 함께M_SQRT1_2= √0,5.

나는 그것을 통계적인 계산에 사용하고 그것은 매우 잘 작동한다.계수를 사용할 필요가 없습니다.

다음은 코드 14줄의 누적 정규 분포의 독립형 C++ 구현입니다.

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

나는 나보다 먼저 대답한 사람들의 제안에 따라 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;
}

여기에 제시된 일반 CDF의 구현은 다음과 같은 단일 정밀도 근사치입니다.float로 대체하다double따라서 7 또는 8개의 유효(최소) 수치만 정확합니다.
Hart의 2배 정밀도 근사 VB 구현은 누적 정규 함수에 대한 West's Better 근사 그림 2를 참조한다.

편집: West 구현의 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;
}

식을 급수 및 연속 분수 근사치에 대해 보다 친숙한 형태로 재배열했습니다.West의 코드의 마지막 매직넘버는 2µ의 제곱근으로, 저는 acos(0) = ½π identity identity identity를 이용하여 컴파일러의 첫 번째 줄에 있는 것을 지연시켰습니다.
매직 넘버를 세 번 체크해 봤지만 뭔가 잘못 입력했을 가능성은 항상 있어요.오타가 발견되면 댓글 달아주세요!

John Cook이 답변에 사용한 테스트 데이터의 결과는 다음과 같습니다.

 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

매스매티카 결과에 대해 주어진 모든 숫자에 동의한다는 사실에서 약간의 위안을 얻습니다.

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.무단 전재 금지.

Boost는 표준:D 여기 있습니다: 수학/통계학을 Boost합니다.

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 대신 2.0을 사용하면 정수 대신 소수점을 얻는 데 도움이 됩니다.

도움이 됐으면 좋겠다.

언급URL : https://stackoverflow.com/questions/2328258/cumulative-normal-distribution-function-in-c-c