具有极小值的函数积分对数的数值稳定评估
Numerically stable evaluation of log of integral of a function with extremely small values
如果我有一个随机数 Z 定义为 sum of two other random numbers、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，因此即使对数概率具有某个有限值，也不能通过简单地取概率的对数来计算它。比较 0.0 的 norm.pdf(39, 1, 0) 和大约 -761 的 norm.logpdf(39, 1, 0)。显然，Scipy 不会将 logpdf 计算为 log(pdf)——它会找到其他方法——因为否则它会 return -inf，一个较差的响应。同样的，我想为我的问题另辟蹊径
（你可能想知道为什么我关心远离均值的值的对数似然。答案是参数拟合。当对数似然是某个巨大的负数时，拟合算法可以更接近，但什么也做不了当它是 -inf 或 nan.)
问题是：有谁知道如何重新排列 log(quad(...)) 以便不计算 quad(...) 从而避免在日志中创建 0.0？
问题是您正在积分的函数的值太小而无法用双精度表示，双精度只能在 1e-308 左右之前使用。
mpmath 助你一臂之力
当双精度不足以进行数值计算时，mpmath, a library for arbitrary-precision floating point operations, is called for. It has its own quad routine, but you'll need to implement your pdf functions so they work at mpmath level (otherwise there won't be anything to integrate). There are many built-in functions, including normal 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，您至少可以尝试 "normalize" 通过将函数的值移动到双精度范围内来尝试 "normalize" 该函数。这是一个例子：
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。
- 正常整合；取对数，对数加上偏移量。从代数上讲，结果只是我们想要取的积分的对数。
如果我有一个随机数 Z 定义为 sum of two other random numbers、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，因此即使对数概率具有某个有限值，也不能通过简单地取概率的对数来计算它。比较 0.0 的 norm.pdf(39, 1, 0) 和大约 -761 的 norm.logpdf(39, 1, 0)。显然，Scipy 不会将 logpdf 计算为 log(pdf)——它会找到其他方法——因为否则它会 return -inf，一个较差的响应。同样的，我想为我的问题另辟蹊径
（你可能想知道为什么我关心远离均值的值的对数似然。答案是参数拟合。当对数似然是某个巨大的负数时，拟合算法可以更接近，但什么也做不了当它是 -inf 或 nan.)
问题是：有谁知道如何重新排列 log(quad(...)) 以便不计算 quad(...) 从而避免在日志中创建 0.0？
问题是您正在积分的函数的值太小而无法用双精度表示，双精度只能在 1e-308 左右之前使用。
mpmath 助你一臂之力
当双精度不足以进行数值计算时，mpmath, a library for arbitrary-precision floating point operations, is called for. It has its own quad routine, but you'll need to implement your pdf functions so they work at mpmath level (otherwise there won't be anything to integrate). There are many built-in functions, including normal 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，您至少可以尝试 "normalize" 通过将函数的值移动到双精度范围内来尝试 "normalize" 该函数。这是一个例子：
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。
- 正常整合；取对数，对数加上偏移量。从代数上讲，结果只是我们想要取的积分的对数。