如果我有一个随机数Z定义为另外两个随机数X和Y的总和，那么Z的概率分布是X和Y的概率分布的卷积。卷积基本上是分布函数的乘积的积分。卷积中的积分通常没有解析解，因此必须使用基本的正交算法进行计算。在伪代码中：
prob_z(z) = integrate(lambda t: prob_x(t) * prob_y(z-t), -inf, inf)
举个具体的例子，正态分布变量X和对数正态分布变量Y的总Z可以用以下 Python/Scipy 代码计算：
from scipy.integrate import quad
from scipy.stats import norm, lognorm
from scipy import log
prob_x = lambda x: norm.pdf(x, 0, 1) # N(mu=0, sigma=1)
prob_y = lambda y: lognorm.pdf(y, 0.1, scale=10) # LogN(mu=log(10), sigma=0.1)
def prob_z(z):
return quad(lambda t: prob_x(t)*prob_y(z-t), -inf, inf)
现在我想计算对数概率。天真的解决方案是简单的做：
def log_prob_z(z):
return log(prob_z(z))
但是，这在数值上是不稳定的。在大约 39 个标准差之后，概率分布在数值上为 0.0，因此即使对数概率具有一些有限值，也不能通过简单地取概率的对数来计算超出此值。比较norm.pdf(39, 1, 0)的 0.0 和norm.logpdf(39, 1, 0)大约 -761。显然，Scipy不会像log(pdf)那样计算logpdf- 它找到其他方式 - 因为否则它会返回-inf，一个较差的响应。同样，我想为我的问题找到另一种方法。
(你可能想知道为什么我关心值的对数相似性，与平均值相去甚远。答案是参数拟合。当对数似然是某个巨大的负数时，拟合算法可以更接近，但当它-inf或nan时，什么也做不了。
问题是：有谁知道如何重新排列log(quad(...))这样我就不会计算quad(...)从而避免在日志中创建 0.0？
问题是你正在集成的函数的值太小，无法以双精度表示，这只在 1e-308 左右之前才有效。
MPMATH救援
当双精度不足以进行数值计算时，需要 mpmath，一个用于任意精度浮点运算的库。 它有自己的quad例程，但你需要实现你的pdf函数，以便它们在mpmath级别工作(否则不会有任何要集成的东西)。有许多内置函数，包括普通的pdf，所以我将用它来说明。
在这里，我使用 SciPy 在相距 70 的距离处卷积两个普通 pdf：
z = 70
p = quad(lambda t: norm.pdf(t, 0, 1)*norm.pdf(z-t, 0, 1), -np.inf, np.inf)[0]
可悲的是，p正好是0.0。
在这里，我在import mpmath as mp之后对 mpmath 做同样的事情：
z = 70
p = mp.quad(lambda t: mp.npdf(t, 0, 1)*mp.npdf(z-t, 0, 1), [-mp.inf, mp.inf])
现在 p 是一个 mpmath 对象，打印为 2.95304756048889e-543，远远超出了双精度刻度。它的对数mp.log(p)是-1249.22086778731。
基于 SciPy 的替代方案：对数偏移
如果由于某种原因无法使用 mpmath，则至少可以尝试通过将函数的值移动到双精度范围内来"规范化"函数。下面是一个示例：
z = 70
offset = 2*norm.logpdf(z/2, 0, 1)
logp = offset + np.log(quad(lambda t: np.exp(norm.logpdf(t, 0, 1) + norm.logpdf(z-t, 0, 1) - offset), -np.inf, np.inf)[0])
这里logp打印-1264.66566393，它不如mpmath结果好(所以我们失去了一些功能)，但它是合理的。我所做的是：
- 计算函数对数最大值的对数(这是变量偏移量)
- 从 PDF 的对数中减去此偏移量;这是
norm.logpdf(t, 0, 1) + norm.logpdf(z-t, 0, 1) - offset - 将结果求幂，因为我们不能只将对数放在积分中。代数上，这与pdfs乘以exp(-offset)的乘积相同。但从数字上看，这是一个不太可能溢出的数字;实际上，在 T = z/2 时，它是 exp(0)=1。
- 正常积分;取对数，对数加偏移量。代数上，结果只是我们想要取的积分的对数。