Calculating the Cumulative Normal Distribution

by M.A. (Thijs) van den Berg


The cumulative normal distribution is defined by the integral

N(x) = \int_{-\infty}^x n(t) dt

with n(t) the normal density function

n(x) = \frac{1}{\sqrt{2\pi}}e^{ -\frac{1}{2} x^2}


The cumulative normal distribution function.
The cumulative normal distribution function.


This integral has no closed form solution, but many good approximations exist.

Contents

[edit] Abromowitz and Stegun approximation

The most used algorithm is algorithm 26.2.17 from Abromowitz and Stegun, Handbook of Mathematical Functions. It has a maximum absolute error of 7.5e^-8.

\begin{align} N(x) &= 1-n(x) \left( b_1t + b_2t^2 + b_3t^3 + b_4t^4 + b_5t^5 \right) + \epsilon \\ t &= \frac{1}{1 + 0.2316419 x}\\ \end{align}

with

b_1 = 0 .31938 1530
b_2 = -0 .35656 3782
b_3 = 1 .78147 7937
b_4 = -1 .82125 5978
b_5 = 1 .33027 4429

[edit] C++ code 

double N(const double x)
{
  const double b1 =  0.319381530;
  const double b2 = -0.356563782;
  const double b3 =  1.781477937;
  const double b4 = -1.821255978;
  const double b5 =  1.330274429;
  const double p  =  0.2316419;
  const double c  =  0.39894228;

  if(x >= 0.0) {
      double t = 1.0 / ( 1.0 + p * x );
      return (1.0 - c * exp( -x * x / 2.0 ) * t *
      ( t *( t * ( t * ( t * b5 + b4 ) + b3 ) + b2 ) + b1 ));
  }
  else {
      double t = 1.0 / ( 1.0 - p * x );
      return ( c * exp( -x * x / 2.0 ) * t *
      ( t *( t * ( t * ( t * b5 + b4 ) + b3 ) + b2 ) + b1 ));
    }
}

[edit] VBA (Visual Basic for Applications) code 

Public Function NormProb(x As Double) As Double

  Dim t As Double
  Const b1 = 0.31938153
  Const b2 = -0.356563782
  Const b3 = 1.781477937
  Const b4 = -1.821255978
  Const b5 = 1.330274429
  Const p = 0.2316419
  Const c = 0.39894228

  If x >= 0 Then
      t = 1# / (1# + p * x)
      NormProb = (1# - c * Exp(-x * x / 2#) * t * _
      (t * (t * (t * (t * b5 + b4) + b3) + b2) + b1))
  Else
      t = 1# / (1# - p * x)
      NormProb = (c * Exp(-x * x / 2#) * t * _
      (t * (t * (t * (t * b5 + b4) + b3) + b2) + b1))
  End If
End Function

[edit] High Precision Approximation

A high (double) precision approximation is described in the article Better Approximations to Cumulative Normal Fuctiona by Graeme West. The algorithm is based on Hart's algorithm 5666 (1968).

Retrieved from "http://www.sitmo.com/doc/Calculating_the_Cumulative_Normal_Distribution"