ECE 150 Fundamentals of programming

The standard normal distribution is defined by the probability density function

$\frac{1}{\sqrt{2 \pi}} e^{-\frac{x^2}{2}}$.

This function does not have a simple anti-derivative, so if we wish to compute the integral, we must use the error function ($erf(x)$):

$\int_{-\infty}^x \frac{1}{\sqrt{2 \pi}} e^{-\frac{\xi^2}{2}} d \xi = \frac{1}{2}\left(1 + erf\left(\frac{x}{\sqrt{2}}\right)\right)$.

Fortunately, this function is implemented in the standard mathematics library:

double std_normal_pdf( double x );
double std_normal_cdf( double x );

double std_normal_pdf( double x ) {
	double const PI{std::acos( -1 )};
	return std::exp( -0.5*x*x )/std::sqrt(2*PI);
}

double std_normal_cdf( double x ) {
	return 0.5 + 0.5*std::erf( x/std::sqrt(2.0) );
}

Now, we often must find when the integral of the standard normal equals a specified value. Unfortunately, this has no solution, so we must approximate the solution using Newton's method:

double inv_int_std_normal( double x );

double std_normal_icdf( double x ) {
	assert( (0 < x) && (x < 1) );
	double result{x - 0.5};

	double s{x - 0.5};
	double ss{s*s};

	double x0;
	double x1{s*(2.50662827463100051
               + ss*(2.62493499095373663
               + ss*(5.77253353861173460
               + ss*(15.6676089632851799
               + ss*47.0357874801131660))))};

	do {
		x0 = x1;
		x1 = x1 - std_normal_cdf( x1 )/std_normal_pdf( x1 );
	} while ( std::abs( x0 - x1 ) > 1e-15 );

	return x1;
}