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).

[edit] Comments

Name (required):

Website:

Comment:


[edit] Amborag said ...

There is a bracket missing in the C++ source code. After the "else" and before the "double". Other than that nice code

--Amborag 13:38, 8 May 2007 (CEST)

[edit] Admin said ...

Thanks!

--Admin 23:52, 20 May 2007 (CEST)

[edit] GuildMaster said ...

This code seems to be for a specific mean and standard deviation. How would I alter it to take in these as parameters?

GMA

--GuildMaster 18:16, 5 June 2007 (CEST)

[edit] GuildMaster said ...

After further thought, if I have mu and sigma for my normal distribution, am I correct that calling this with (x/sigma - mu/sigma) would accomplish what I want?

GMA

--GuildMaster 20:32, 5 June 2007 (CEST)

[edit] Admin said ...

Yes this function is for zero mean and a standard deviation of 1.

Your solution for arbitrary mean (mu) and standard deviation (sigma) is correct.

--Admin 16:21, 26 June 2007 (CEST)

[edit] L said ...

Concerning GuildMaster's question: Are you sure? Shouldn’t all expression be affected by “*1/sigma” too? L

--L 12:52, 29 June 2007 (CEST)

[edit] Afi said ...

I have a problem with number other than(mean=0 and sd=1). Please say me what can i do till i can find the number of norm as excel calculate.

--Afi 21:37, 27 August 2007 (CEST)

[edit] sitmo said ...

N(x,mean,sigma) = N( (x-mean)/sigma )

--sitmo 12:29, 6 September 2007 (CEST)

[edit] Michael Andrews said ...

The link you have provided also provides an approximation to n(x) (in the section linked to it's called Z(x)). By computing with the exp function I think your method might not have the same upper bound on error (\epsilon).

--Michael Andrews 22:40, 2 October 2007 (CEST)