当'a'和'b'远离均值时，Python和R中的截断正态分布

Question

当尝试获取截断正态分布的概率密度函数（pdf）时：

from scipy.stats import truncnorm
truncnorm.pdf(-31, np.inf, -30, loc=0, scale=1)

这个过程运行得很好。但是如果上限离均值太远，分配给非截断部分（整体质量应该是1的地方）的样本概率会变成NaN：

# -41 is one of the points with highest probability. Why nan?
>truncnorm.pdf(-41, np.inf, -40, loc=0, scale=1)
nan

# 39 is impossible since it lays in the truncated side
>truncnorm.pdf(-39, np.inf, -40, loc=0, scale=1)
0.0

这是因为数值精度的问题还是其他什么原因导致的bug吗？有没有其他方法可以解决这个问题？

更新 1（使用R库“truncnorm”）：

这似乎是一个常见的问题。在R的“truncnorm”库中也遇到了同样的问题：

 > dtruncnorm(-41, a=-Inf, b=-40, mean = 0, sd = 1)
 [1] NaN

更新 2（使用R库“msm”）：

在他的博客中，Christian Robert 指出了“msm”库，它实现了他的论文。

然而，在这种情况下它也出现了崩溃：

> dtnorm(-41, mean = 0, sd=1, lower=-Inf, upper=-40)
[1] NaN

Answer 1

truncnorm的计算是基于正态分布的累积分布函数（cdf）。

目前为止，在分布的尾部，无法用浮点数（双精度）精确表示这个累积分布函数。

>>> stats.norm.cdf(-37)
5.7255712225239266e-300
>>> stats.norm.cdf(-38)
0.0

>>> stats.norm.pdf(-37)
2.120006551524606e-298
>>> stats.norm.pdf(-38)
1.0972210519949712e-314
>>> stats.norm.pdf(-39)
0.0

>>> np.finfo(float).tiny
2.2250738585072014e-308

实现这个的唯一方法就是直接计算或近似这个截断分布，而不是通过正态分布的特殊函数。

我从来没有见过需要使用这个的情况。