# Moments of the truncated normal distribution (univariate) away from the mean

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:

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

• @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. – becko Dec 5 '17 at 12:08
• 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. – whuber Dec 6 '17 at 18:32
• 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. – whuber Dec 6 '17 at 19:29
• 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. – whuber Dec 6 '17 at 20:35
• @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]. – whuber Jul 15 '19 at 20:04

## 1 Answer

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?