Direct computation of $\operatorname{log}(\operatorname{cdf})$ for a normal distribution

This question is linked to the normal distribution for a random variable. The probability density function (pdf) is expressed as:

$$\frac{1}{\sigma\sqrt{2\pi}}\, e^{-\frac{(x - \mu)^2}{2 \sigma^2}}$$

and the cumulative distribution function (cdf) is expressed as:

$$\frac12\left[1 + \operatorname{erf}\left( \frac{x-\mu}{\sigma\sqrt{2}}\right)\right]$$

pdf or cdf functions might take extremely low values. When computing logarithm of pdf,

$$- \operatorname{log}(\sigma\sqrt{2\pi}) - \frac{(x - \mu)^2}{2 \sigma^2}$$

is used directly, instead of first computing eq. 1, then taking the log.

Currently I compute the log of the cdf by first computing eq. 2, then taking its log. This creates numerical underflows in my program. I would like to compute the log directly, using a formula similar to eq. 3.

Is it possible ?

• I think this depends on what your numerical problem is: is it full underflows (numbers going below $2^{-1024}$) or is it roundoff errors (numbers on the order of $1$ getting added to numbers on the order of $2^{-53}$)? The solutions to these two kinds of problems are rather different. – Ian Jun 30 '15 at 15:45
• It is a bit of both. I need to add various log. The majority of which are close to 0., but some are under 1.7E- 308, which is the limit of the double type I use to compute the cdf. When I compute the sum I get $-\infty$. – vkubicki Jun 30 '15 at 15:51
• When adding numbers of enormously varying magnitude, it's generally a good idea to find the summand of largest magnitude, divide all the terms by that, add them, and then multiply by what you divided by. – Ian Jun 30 '15 at 16:30
• I do this in other part of the code. The problem here is that can not even compute the value of log(cdf). – vkubicki Jun 30 '15 at 17:20
• You can't compute individual summands? How small should your actual answer be, more or less? It seems that if you can't compute individual summands and all of your summands have the same sign then your answer might not even be a floating point representable number. In this case you would need to use a higher precision arithmetic library to do anything. – Ian Jun 30 '15 at 17:23