4

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

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

举一个具体的例子，Z一个正态分布变量X和一个对数正态分布变量之和Y可以用下面的 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。logpdf显然，Scipy 不会计算log(pdf)- 它会找到其他方式 - 因为否则它会返回-inf，一个较差的响应。同样，我想为我的问题找到另一种方法。

（您可能想知道为什么我关心与平均值相差甚远的值的对数相似度。答案是参数拟合。当对数似然为某个极大的负数时，拟合算法可以更接近，但当它是-inf或时，什么都做不了nan。 )

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

4

1 回答 1

4

问题是您正在积分的函数的值太小而无法以双精度表示，这仅在 1e-308 左右之前是好的。

mpmath 来救援

当双精度不足以进行数值计算时，需要一个用于任意精度浮点运算的库mpmath 。它有自己的quad例程，但您需要实现您的 pdf 函数，以便它们在 mpmath 级别工作（否则将没有任何东西可以集成）。有很多内置函数，包括普通的 pdf，所以我将使用它来进行说明。

在这里，我将两个相距 70 的普通 pdf 与 SciPy 进行卷积：

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。
  • 正常整合；取对数，对数加上偏移量。在代数上，结果只是我们想要取的积分的对数。
于 2017-11-30T02:57:57.387 回答