scipy wired output of scipy.stats.truncnorm

Question

I am using scipy version of random generation for truncated normal distribution:

from scipy.stats import truncnorm; truncnorm.rvs(8.26,1e10)

if I try multiple times, the supposed random generation would freeze to 1.7976931348623157e+308, which is wrong...

Can anyone help me to understand what happened here?

Thanks

Answer 1

truncnorm uses the ppf (inverse function of the cdf) to transform the random variable.

It looks like the way the ppf is calculated doesn't work when the range is so far out in the tail

>>> truncnorm._ppf(np.linspace(0, 1, 11), 8.26,1e10)
array([ 8.20953615,  8.20953615,         inf,  8.20953615,  8.20953615,
               inf,         inf,         inf,         inf,         inf,
               inf])

The problem is that both values for the truncation are far out in the tail of the normal distribution, so that we are essentially only working with 1 minus rounding problems. That is a floating point restriction, that we don't have higher precision than double, for example 1-1e-20 == 1.

>>> stats.norm.cdf(8.26)
0.99999999999999989

However, the normal distribution is symmetric and at the lower tail the ppf doesn't have the floating point problems because we are working in the neighborhood of zero. I think the following should give correct results

>>> -truncnorm.rvs(-1e10, -8.26,size=10)
array([ 8.8359488 ,  8.3112093 ,  8.36978251,  8.5732518 ,  8.57420297,
        8.27919164,  8.5692663 ,  8.28482223,  8.35149422,  8.47994703])