Mathematica使用符号积分（不是近似值！），报告一个等于 1.6534640367102553437 到 20 位十进制数字的值。

red = GammaDistribution[20/17, 17/20];
gray = InverseGaussianDistribution[1, 832/1000];
kl[pF_, qF_] := Module[{p, q},
   p[x_] := PDF[pF, x];
   q[x_] := PDF[qF, x];
   Integrate[p[x] (Log[p[x]] - Log[q[x]]), {x, 0, \[Infinity]}]
   ];
kl[red, gray]

一般来说，使用小蒙特卡罗样本不足以计算这些积分。正如我们在另一个线程中看到的那样，KL 散度的值可以由基础 PDF 非零且另一个 PDF 接近于零的短区间内的积分决定。小样本可能会完全错过这样的间隔，并且可能需要大样本才能达到足够的时间才能获得准确的结果。

看一下 0 附近逆高斯的 Gamma PDF（红色，虚线）和对数 PDF（灰色）：

在这种情况下，伽玛 PDF 保持在 0 附近的高位，而逆高斯 PDF 的对数在那里发散。几乎所有的 KL 散度都是由区间 [0, 0.05] 中的值贡献的，在 Gamma 分布下该值的概率为 3.2%。 Gamma 变量的样本中落入此区间的元素数量因此具有二项式 (.032, ) 分布。它的标准差等于 = 0.55%。因此，作为粗略估计，我们不能期望您的积分具有比 (3.2% - 2*0.55%) / 3.2% = 大约 40% 更好的相对准确度，因为它采样的次数关键间隔有这么多的错误。这解释了你的结果和$N = 1000$$N$$.032(1 - .032)/\sqrt{N}$数学的. 要将此误差降低到 1%（仅精确到小数点后两位），您需要将样本大约乘以：也就是说，您需要几百万个值。$(40/1)^2 = 1600$

这是KL 散度的 1000 个独立估计的自然对数的直方图。每个估计值平均从 Gamma 分布中随机获得 1000 个值。正确的值显示为红色虚线。

Histogram[Log[Table[
 Mean[Log[PDF[red, #]/PDF[gray, #]] & /@ RandomReal[red, 1000]], {i,1,1000}]]]

虽然平均而言，这种蒙特卡洛方法是无偏的，但大多数时候（在这个模拟中为 87%）估计值太低了。为了弥补这一点，有时高估可能是粗略的：这些估计中最大的是 18.98。（广泛分布的值表明 1.286916 的估计实际上没有可靠的十进制数字！）由于分布中的这种巨大偏斜，情况实际上比我之前用简单的二项式思维实验估计的要糟糕得多。这些模拟的平均值（总共包含 1000*1000 个值）仅为 1.21，仍比真实值低约 25%。

一般来说，为了计算 KL 散度，您需要使用自适应正交或精确方法。