11
$\begingroup$

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.

$\endgroup$
5
  • $\begingroup$ How far below $0$ could $x$ go? Different useful approximations kick in below around $-4$ to $-5$ and yet other ones are needed for very negative values. For instance, if $x\lt -38$ or so, the CDF will underflow IEEE double precision floats (but of course its log will not be problematic). How accurate does your approximation need to be? $\endgroup$
    – whuber
    Jul 6, 2014 at 23:15
  • $\begingroup$ @whuber x won't go below -5. Would be great if the log_2 could be within 0.2 + 10% or so of its true value, but this isn't critical. I'd probably go with whatever I can get. $\endgroup$
    – Museful
    Jul 6, 2014 at 23:45
  • 1
    $\begingroup$ That's rather critical information in your comment there. It's probably better to also add it to your post. $\endgroup$
    – Glen_b
    Jul 7, 2014 at 0:19
  • $\begingroup$ Related: stats.stackexchange.com/q/7200/2970. The Cody and/or Lee references may provide suitable approaches. $\endgroup$
    – cardinal
    Jul 7, 2014 at 3:19
  • $\begingroup$ @Cardinal That's right: the Mills ratio yields excellent approximations to the logarithm for extremely negative $x$. By empirically adjusting its coefficients one can achieve accuracy of $\pm 0.0003$ (in the natural logs) for $x\in [-7,-3]$. By taking the next term in the continued fraction and again adjusting the coefficients, errors can be limited to $\pm 0.000005$ for $x\lt -7$: that approximation is $\log(R(x))\approx-\log \left(x+\frac{1}{x+\frac{1.9}{x}}\right)$. But the more terms of the continued fraction you take, the larger $|x|$ must be for it to be a good approximation. $\endgroup$
    – whuber
    Jul 7, 2014 at 15:10

8 Answers 8

8
$\begingroup$

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:

enter image description here

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

enter image description here

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.

enter image description here

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

enter image description here

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.

$\endgroup$
5
  • 1
    $\begingroup$ For a number of reasons, that Wikipedia page needs some significant attention, or, at least, some TLC from a knowledgeable statistician or probabilist. While not the only example, the statement regarding the K&L approximation would make any mathematician cringe since the absolute relative error is obviously (i.e., by inspection) nearly 90% in the tails! :-) $\endgroup$
    – cardinal
    Jul 7, 2014 at 3:01
  • $\begingroup$ (+1) Simple (e.g., polynomial or rational) approximations can do quite well over a wide range of values of the argument for $\log \Phi$ since it's actually quite linear near zero. A simple graph would depict this nicely. $\endgroup$
    – cardinal
    Jul 7, 2014 at 3:04
  • $\begingroup$ @cardinal Thanks for both comments. I've added some graphs, including one for $\log_2(\Phi(x))$ (in which the near-linearity at zero can be seen), another for the error in the quintic approximation I got, and I also had a closer look at the K&L approximation on the log-scale (again with graphs). $\endgroup$
    – Glen_b
    Jul 7, 2014 at 4:28
  • $\begingroup$ I didn't expect such a comprehensive answer! By "<0.2+10% or so" I meant $<0.2+0.1*\log_2(\Phi(x))$ i.e. better accuracy closer to $x=0$, but by pre-constraining the polynomial's y-intercept to the correct value I'm getting exactly what I needed. Thank you so much! $\endgroup$
    – Museful
    Jul 7, 2014 at 9:44
  • $\begingroup$ I'm glad it was useful to you. $\endgroup$
    – Glen_b
    Jul 7, 2014 at 10:03
4
$\begingroup$

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))
$\endgroup$
3
  • 1
    $\begingroup$ Cool trick. Shouldn't this be -F.softplus(np.log(2) - np.sqrt(2*np.pi)*x) though? $\endgroup$ Jun 23, 2020 at 19:16
  • $\begingroup$ @huehue Yes! I just corrected the formula, well caught. $\endgroup$ Jun 24, 2020 at 11:53
  • 1
    $\begingroup$ If anyone wants a really accurate and numerically stable PyTorch implementation for log_norm_cdf (for the entire range), see my gist: gist.github.com/chausies/011df759f167b17b5278264454fff379 $\endgroup$
    – chausies
    Jun 30, 2021 at 11:29
2
$\begingroup$

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

$\endgroup$
2
$\begingroup$

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.

$\endgroup$
1
  • $\begingroup$ This doesn't seem to be of any help, because it starts with a computation that is even more complicated than the function to be approximated! $\endgroup$
    – whuber
    Jun 30, 2017 at 15:17
2
$\begingroup$

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

enter image description here

$\endgroup$
1
  • 1
    $\begingroup$ +1. This is a good approach for $x \ll 0$. It's called the "Mills ratio". Please note that comments to the question indicate there is originally no interest in such extreme values: it focuses on $-5\le x \lt 0$. Your plots disguise considerable error for smaller values of $|x|$. $\endgroup$
    – whuber
    Jun 30, 2017 at 15:19
0
$\begingroup$

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

$\endgroup$
0
$\begingroup$

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

enter image description here

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.

$\endgroup$
0
$\begingroup$

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

```
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.