具有极小值的函数积分对数的数值稳定评估

Question

如果我有一个随机数Z，它被定义为另外两个随机数 和 的总和，X那么Y的概率分布是和Z的概率分布的卷积。卷积基本上是分布函数乘积的积分。卷积中的积分通常没有解析解，因此必须使用基本求积算法来计算。在伪代码中：XY

\n\n

prob_z(z) = integrate(lambda t: prob_x(t) * prob_y(z-t), -inf, inf)\n

\n\n

举一个具体的例子，可以使用以下Python/Scipy代码计算Z正态分布变量X和对数正态分布变量的总和：Y

\n\n

from scipy.integrate import quad\nfrom scipy.stats import norm, lognorm\nfrom scipy import log\n\nprob_x = lambda x: norm.pdf(x, 0, 1)  # N(mu=0, sigma=1)\nprob_y = lambda y: lognorm.pdf(y, 0.1, scale=10)  # LogN(mu=log(10), sigma=0.1)\ndef prob_z(z):\n    return quad(lambda t: prob_x(t)*prob_y(z-t), -inf, inf)\n

\n\n

现在我想计算对数概率。天真的解决方案是简单地执行以下操作：

\n\n

def log_prob_z(z):\n    return log(prob_z(z))\n

\n\n

然而，这在数值上是不稳定的。在大约 39 个标准差之后，概率分布在数值上为 0.0，因此即使对数概率具有某个有限值，也不能通过简单地取概率的对数来计算此之外的值。比较norm.pdf(39, 1, 0)哪个是 0.0norm.logpdf(39, 1, 0)哪个是 -761。显然，Scipy 不会计算logpdf为log(pdf)\xe2\x80\x94，它会找到其他方式 \xe2\x80\x94，因为否则它会返回-inf，这是一个较差的响应。同样，我想为我的问题找到另一种方法。

\n\n

（你可能想知道为什么我关心远离平均值的值的对数似然。答案是参数拟合。当对数似然是某个巨大的负数时，拟合算法可以变得更接近，但当它是 或 时，什么也做不-inf了nan。 ）

\n\n

问题是：有谁知道我如何重新排列，log(quad(...))这样我就不会计算quad(...)，从而避免在日志中创建 0.0 ？

\n

Answer 1

问题是您要积分的函数值太小，无法以双精度表示，而双精度仅在 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。

在这里我对 mpmath 做了同样的事情，之后import mpmath as mp：

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。
• 正常整合；取对数，在对数上加上偏移量。从代数上来说，结果就是我们想要取的积分的对数。