Approximation of logarithm of standard normal CDF for x<0

Question

Does anyone know of an approximation for the logarithm of the standard normal CDF for x<0?

I need to implement an algorithm that very quickly calculates it. The straightforward way, of course, is to first calculate the CDF (for which I can find suitably simple approximations on Wikipedia), and then to take the logarithm of that. Obviously I'd like to avoid the time-cost of having to calculate two special functions, not to mention that the tiny intermediate value (in the tail) rules out using fixed-point arithmetic which is much faster than floating-point arithmetic.

I know there are hundreds of approximations for all kinds of statistical functions, but the fact that this is the logarithm of one makes it harder to find one. I'd be very grateful if anyone could point me to one, or to a source where I might find one.

Answer 1

Even a simple least squares cubic fit to $\log(\Phi(x))$ values for $x$ between -5 and 0 seems to be fairly adequate.

Just to get a quick sense of it I generated $x$ values every 0.01 between -5 and 0, and tried least squares cubic and quintic (5th degree) polynomial fits to $\log_2(\Phi(x))$. I presume you can do that as easily as I did, so I won't labor the point.

The maximum absolute error in $\log_2(\Phi(x))$ for the quintic is $5.2\times 10^{-4}$, which occurs at zero.

[It's not completely clear what you mean by "within 0.2 + 10% or so". If you could elaborate so as to make your criterion explicit, then I could address that in detail, and maybe adjust the weights to better optimize your criterion.]

When evaluating polynomials, if speed matters, you should keep Horner's method in mind.

---

As cardinal suggested, the plot is quite illustrative. Since I'd already generated one, I should have put it here:

Here's the (signed) error (absolute scale error in the logs):

it looks much as we might expect for a least squares fit. The lack of fit is pretty moderate; for many purposes that would be fine.

---

An alternative is you can flip the Karagiannidis & Lioumpas approximation (see here) for the upper tail around to the lower tail (by replacing $x$ by $-x$ in their formula) and taking logs:

$$Q(x)\approx\frac{\left( 1-e^{-1.4x}\right) e^{-\frac{x^{2}}{2}}}{1.135\sqrt{2\pi}x}, x >0 $$

So we get

$$\ln(\Phi(x))\approx \ln(1-e^{1.4x})-\ln(-x) -\frac{x^2}{2} - 1.04557$$

to get base 2 logs, you just multiply the result by $\log_2(e)$

This is less accurate than the quintic fit I mentioned, and likely slower to evaluate because of the logs and exponentiation. Still, it's nice and short, which has some advantages. Note that it wasn't designed for the log-scale.

---

The original paper has K&L's Equation 6 as:

$$\frac{(1-e^{-Ax}) e^{-x^2}}{B\sqrt{\pi}x}\approx\text{erfc}(x)$$

For $x$ values on the range [0,20] they suggest $A=1.98$ and $B=1.135$.

Dividing $x$ through by $\sqrt{2}$ to obtain $Q$, that suggests the formula on the Wikipedia page ($1.98/\sqrt{2}\approx 1.40007$)

---

Let's look at the quality of the approximation on the log-scale.

The purple dashed line is the K&L approximation. Now let's look at the error:

The size of the error gets quite a bit larger than our quintic. The fact that it's nearly linear in $x$ at the left half of the range suggests that we might improve the error by adding about $0.02x$ or $0.025x$ to the formula in the logs - but this would make the approximation a little worse for values near -1. Whether one would do that depends on what characteristics are desired.

Answer 2

Here is an alternate solution for neural network users:

You can approximate the CDF with 1 / (1 + 2*exp(-sqrt(2*pi)*x)) (extending this approximation).

And exploit the fact that most deep learning framework have a numerically stable implementation of log(1 + exp(x)) as the softplus activation function.

The resulting formula can be written in two lines of pytorch:

def log_standard_normal_cdf(x):
    return -F.softplus(np.log(2) - x*np.sqrt(2*np.pi))

Answer 3

To cover a wider range than originally requested, I eventually came up with this rational approximation

$$ \ln(\Phi(x<0)) \approx -\frac{1}{2}x^2 -4.8 + 2509\frac{x-13}{(x-40)^2(x-5)} $$

which has absolute error under 0.04 out to 20 standard deviations, after which the error may exceed that but remains under 0.04% of the value (as far as I can verify with numerical calculations).

Answer 4

I recently had the same problem and found what I think is a nicer solution. The Faddeeva package has a set of accurate functions for computing erf(z), erfc(z), etc. They all start by computing the Faddeeva function w(z):

\begin{equation} w(z) = e^{-z^2} \text{erfc}(-iz) \end{equation}

you can show that the log of the normal cdf is

\begin{eqnarray} \log\Phi(z) &= \log\left[ \frac{1}{2}\left( 1 + \text{erf}\left( \frac{z}{\sqrt{2}} \right) \right) \right] \\ &= \log\left(\frac{1}{2}\right) + \log\left[ 1 + \text{erf}\left(\frac{z}{\sqrt{2}} \right) \right] \\ &= \log\left(\frac{1}{2}\right) + \log\left[e^{-z^2} w\left(\frac{-iz}{\sqrt{2}}\right) \right] \\ &= \log\left(\frac{1}{2}\right) - \frac{z^2}{2} + \log\left[ w\left(\frac{-iz}{\sqrt{2}}\right) \right] \end{eqnarray}

So now we can compute three easy quantities that are all nicely behaved and get $\log\Phi(z)$ also nicely behaved for very small z. I verified that you get exactly the same values as a regular normcdf for the domain in which it is valid.

For very small z, I don't know the "right" answer, but for z=-100 you get -5005.5242 which normally would just give you -inf. And this method has the advantage of using existing packages for doing the tricky computations so you don't have to make up your own approximations and it is fast.

In python, $w(z)$ is available as scipy.special.wofz.

Answer 5

We can use l'Hôpital's Rule:

$$ \lim_{x \rightarrow -\infty} \frac{\Phi(x)}{\phi(x)} = \frac{\phi(x)}{-x\phi(x)}=\frac{1}{-x} $$

So:

$$ \lim_{x \rightarrow -\infty} \Phi(x) = \frac{\phi(x)}{-x} $$

And then we get:

$$ \lim_{x \rightarrow -\infty} \log \Phi(x) = \log(\frac{1}{\sqrt{2\pi}}) -0.5 x^2 - \log(-x) $$

This approximation should be accurate if x is small.

Also try this matlab code, it gives really good approximation: I found it in the GPML toolbox

% Safe computation of logphi(z) = log(normcdf(z)) and its derivatives
%                    dlogphi(z) = normpdf(x)/normcdf(x).
% The function is based on index 5725 in Hart et al. and gsl_sf_log_erfc_e.
%
% Copyright (c) by Carl Edward Rasmussen and Hannes Nickisch, 2013-11-13.
function [lp,dlp,d2lp,d3lp] = logphi(z)
  z = real(z);                                 % support for real arguments only
  lp = zeros(size(z));                                         % allocate memory
  id1 = z.*z<0.0492;                                 % first case: close to zero
  lp0 = -z(id1)/sqrt(2*pi);
  c = [ 0.00048204; -0.00142906; 0.0013200243174; 0.0009461589032;
       -0.0045563339802; 0.00556964649138; 0.00125993961762116;
       -0.01621575378835404; 0.02629651521057465; -0.001829764677455021;
       2*(1-pi/3); (4-pi)/3; 1; 1];
  f = 0; for i=1:14, f = lp0.*(c(i)+f); end, lp(id1) = -2*f-log(2);
  id2 = z<-11.3137;                                    % second case: very small
  r = [ 1.2753666447299659525; 5.019049726784267463450;
        6.1602098531096305441; 7.409740605964741794425;
        2.9788656263939928886 ];
  q = [ 2.260528520767326969592;  9.3960340162350541504;
       12.048951927855129036034; 17.081440747466004316; 
        9.608965327192787870698;  3.3690752069827527677 ];
  num = 0.5641895835477550741; for i=1:5, num = -z(id2).*num/sqrt(2) + r(i); end
  den = 1.0;                   for i=1:6, den = -z(id2).*den/sqrt(2) + q(i); end
  e = num./den; lp(id2) = log(e/2) - z(id2).^2/2;
  id3 = ~id2 & ~id1; lp(id3) = log(erfc(-z(id3)/sqrt(2))/2);  % third case: rest
  if nargout>1                                        % compute first derivative
    dlp = zeros(size(z));                                      % allocate memory
    dlp( id2) = abs(den./num) * sqrt(2/pi); % strictly positive first derivative
    dlp(~id2) = exp(-z(~id2).*z(~id2)/2-lp(~id2))/sqrt(2*pi); % safe computation
    if nargout>2                                     % compute second derivative
      d2lp = -dlp.*abs(z+dlp);             % strictly negative second derivative
      if nargout>3                                    % compute third derivative
        d3lp = -d2lp.*abs(z+2*dlp)-dlp;     % strictly positive third derivative
      end
    end
  end

strong text

When x is very small,like x = -100, direct calculation of $\log\Phi(x)$ is impossible,but the approximation is still accurate

Answer 6

Note that it's possible to avoid using expressions such as $\lim_{x\to-\infty}\Phi(x)=-\phi(x)/x$: Using L'Hopital's twice, we see that $$ \lim_{x\to-\infty} -\frac{x\Phi(x)}{\phi(x)} = \lim_{x\to-\infty} \frac{x\phi(x)+\Phi(x)} {x\phi(x)} = 1+ \lim_{x\to-\infty} \frac{\phi(x)}{\phi(x) - x^2\phi(x)} = 1, $$ which means that $$ \Phi(x) \approx -\frac{\phi(x)}{x} $$ is a reasonable approximation for $x\ll0$.

Answer 7

You can also use monotone cubic interpolation using e.g. the Fritsch–Carlson method. This can give you a very precise solution and you can make an implementation like in this package. Here is an example with 300 knots:

# test values
us <- seq(-5, 0, length.out = 10000)

# create the spline function
n_points <- 300L
eps <- 1e-9
x <- seq(qnorm(eps), 1, length.out = n_points)
f <- splinefun(x = x, y = pnorm(x, log.p = TRUE), method = "monoH.FC")

max(abs(pnorm(us, log.p = TRUE) - f(us))) # largest abs error
#R> [1] 4.488093e-08

# plot x versus the error
xs <- seq(-5, 0, length.out = 1000)
plot(xs, pnorm(xs, log.p = TRUE) - f(xs), xlab = "x", ylab = "error", 
     type = "h") 

This R version is not fast. However, you can make a fast implementation like in the package I linked to using that there is a fixed distance between all the knots. This allows you to quickly find the closets knot for a given value.

One though will need to an approximation to plug in for input values smaller than the smallest knot.

Answer 8

I'm kinda late to the party but I would suggest to compute this quantity without approximations, I was using one of the approximations proposed in this post for academic reasons (code for my MSc thesis) and since I really need very high precision I ended up rewriting the function in pure C by using MPFR to avoid numerical problems

#include <gmp.h>
#include <mpfr.h>

double norm_logcdf(double x){
    mpfr_t mpfr_x;
    mpfr_init2(mpfr_x,64);
    mpfr_set_d(mpfr_x,x,MPFR_RNDN);
    mpfr_neg(mpfr_x,mpfr_x,MPFR_RNDN);
    mpfr_mul_d(mpfr_x,mpfr_x, (double) M_SQRT1_2, MPFR_RNDN);
    mpfr_erfc(mpfr_x,mpfr_x,MPFR_RNDN);
    mpfr_mul_d(mpfr_x,mpfr_x,(double) 0.5, MPFR_RNDN);
    mpfr_log(mpfr_x,mpfr_x,MPFR_RNDN);
    double res = mpfr_get_d(mpfr_x,MPFR_RNDN);
    mpfr_clear(mpfr_x);
    return res;
}

```