3
$\begingroup$

I need to compute the mean and variance of the truncated normal distribution. For simplicity, let us focus on a standard normal, since the general case can be reduced to this. The PDF is given by:

$$f\left(x;a,b\right)=\begin{cases} \frac{\phi\left(x\right)}{\Phi\left(b\right)-\Phi\left(a\right)} & a\le x\le b\\ 0 & \text{otherwise} \end{cases}$$

where

$$\phi\left(x\right)=\frac{1}{\sqrt{2\pi}}\mathrm{e}^{-x^{2}/2},\quad\Phi\left(x\right)=\frac{1}{2}\left(1+\mathrm{erf}\left(x/\sqrt{2}\right)\right)$$

The analytical formulas for the mean and variance of $f$ are:

$$\left\langle x\right\rangle =\frac{\phi\left(a\right)-\phi\left(b\right)}{\Phi\left(b\right)-\Phi\left(a\right)}$$

$$\mathrm{var}\,x =1+\frac{a\phi\left(a\right)-\beta\phi\left(b\right)}{\Phi\left(b\right)-\Phi\left(a\right)}-\left(\frac{\phi\left(a\right)-\phi\left(b\right)}{\Phi\left(b\right)-\Phi\left(a\right)}\right)^{2}$$

However, if you try to evaluate those expressions numerically, when $a,b$ are sufficiently away from the mean, you will encounter multiple numerical issues. I have tried Mathematica, Python scipy, Julia Distributions.jl, and they all give NaN.

For example, this is what happens in Mathematica:

enter image description here

Has anyone encountered this issue before? Anyone can suggest a solution, or an alternative package?

Improve this question
$\endgroup$
15
  • $\begingroup$ @COOLSerdash Precisely, it's a numerical issue. As I said, it occurs when $a,b$ are sufficiently away from the mean. For $a=-1, b=1.5$ it works as expected. However I need to evaluate the moments away from the mean. $\endgroup$
    – a06e
    Dec 5, 2017 at 12:08
  • 1
    $\begingroup$ On the contrary, numerical integration will be extremely fast because you can closely approximate the integral with just a few evaluations at arguments close to the minimum of $|a|$ and $|b|$. The reason is that almost all the probability is focused near that minimum. You have no hope of doing these calculations in double precision arithmetic, though. This all implies you should be interested in approximating the moments. How you do so depends on your accuracy needs. $\endgroup$
    – whuber ♦
    Dec 6, 2017 at 18:32
  • 2
    $\begingroup$ A software package is no answer. But you can do better. E.g., the logarithm of the standard Normal tail probability for large $a$ can be very closely approximated as $$\log \Pr(Z \gt a) \approx -\frac{\log(2\pi)}{2}-\frac{a^2}{2}-\log(a)-\frac{1}{a^2}+\frac{5}{2a^4}-\frac{37}{3a^6}.$$ Even for $a$ as small as $4$ this is accurate to three decimal places (that is, $0.1\%$ relative accuracy). You can obtain analogous formulas for the moments. Use linear combinations of these results to obtain the moments of the truncated distribution. $\endgroup$
    – whuber ♦
    Dec 6, 2017 at 19:29
  • 2
    $\begingroup$ Your reference relies on having access to accurate computation of $\operatorname{erf}$ for large arguments--and that's hard to find (and perhaps can be an expensive computation, especially when using extended-precision arithmetic). The Laurent expansion of the logarithm that I gave suffers from none of those problems. For $a \gt 100$ you can't get any more accurate (in double precision arithmetic) and it's hard to imagine a computationally simpler approach. The advantage of using a package is that somebody has already written (and, we hope, tested) it--so use it if it works. $\endgroup$
    – whuber ♦
    Dec 6, 2017 at 20:35
  • 2
    $\begingroup$ @Becko It's a series expansion at infinity. You can get Wolfram Alpha to come close (the result will need some minor tweaking) with Series Log(erfc(x)), (x, infinity, 6) If you're using Mathematica directly, the command is Simplify[ Series[Log[CDF[NormalDistribution[], -a]], {a, Infinity, 6}], Assumptions -> a > 0]. $\endgroup$
    – whuber ♦
    Jul 15, 2019 at 20:04

1 Answer 1

Reset to default
2
$\begingroup$

I found a package in Julia that does this, and seems to take care of truncation away from the mode as far as I have tried it,

https://github.com/cossio/TruncatedNormal.jl

An explanation of why it works is given here. Basically the trick is to substitute the erf by related functions (erfc, erfcx) at different parameter regimes, so as to avoid catastrophic cancellation at the subtractions involved.

Update: Actually the above library has issues computing the variance (though the mean is fine). For example, for a standard normal distribution truncated to the interval $[a,b]$ with $a=-10^6$ and $b=-999000$, this package computes a negative variance:

tnvar(a,b) = -0.0001220703125

which is obviously wrong.

Any ideas?

Improve this answer
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

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