python - 编写集成高斯函数的 Python 函数的最佳方法？

Question

在尝试使用 scipy 的 quad 方法对高斯积分(假设有一种名为 gauss 的高斯方法)时，我在将所需参数传递给高斯并让 quad 对正确变量进行积分时遇到了问题。有没有人有一个很好的例子来说明如何使用带有多维函数的四边形？

但这让我想到了一个更宏大的问题，即整合高斯的最佳方法。我没有在 scipy 中找到高斯积分(令我惊讶)。我的计划是编写一个简单的高斯函数并将其传递给 quad(或者现在可能是一个固定宽度的积分器)。你会怎么办？

编辑:固定宽度的意思是像 trapz 一样使用固定的 dx 来计算曲线下的面积。

到目前为止，我想到的是一个方法 make___gauss，它返回一个 lambda 函数，然后可以进入 quad。通过这种方式，我可以在积分之前使用我需要的平均值和方差制作一个正常函数。

def make_gauss(N, sigma, mu):
    return (lambda x: N/(sigma * (2*numpy.pi)**.5) *
            numpy.e ** (-(x-mu)**2/(2 * sigma**2)))

quad(make_gauss(N=10, sigma=2, mu=0), -inf, inf)

当我尝试传递一个通用的高斯函数(需要用 x、N、mu 和 sigma 调用)并使用 quad 填充一些值时

quad(gen_gauss, -inf, inf, (10,2,0))

参数 10、2 和 0 不一定匹配 N=10、sigma=2、mu=0，这提示了更扩展的定义。

scipy.special 中的 erf(z) 需要我准确定义 t 最初是什么，但很高兴知道它在那里。

Answer 1

好吧，你似乎对几件事很困惑。让我们从头开始:您提到了“多维函数”，然后继续讨论通常的单变量高斯曲线。这不是一个多维函数:当你对它进行积分时，你只对一个变量 (x) 进行积分。区别很重要，因为有一个叫做“多元高斯分布”的怪物，它是一个真正的多维函数，如果集成，需要对两个或多个变量进行积分(使用我之前提到的昂贵的蒙特卡罗技术)。但是您似乎只是在谈论常规的单变量高斯，它更容易使用、集成等等。

一变量高斯分布有两个参数，sigma和 mu ，并且是单个变量的函数，我们将表示为 x .您似乎还携带了标准化参数 n (这在一些应用程序中很有用)。归一化参数通常不包含在计算中，因为您可以在最后重新添加它们(请记住，积分是一个线性运算符:int(n*f(x), x) = n*int(f(x), x))。但如果您愿意，我们可以随身携带；我喜欢的正态分布符号是
N(x | mu, sigma, n) := (n/(sigma*sqrt(2*pi))) * exp((-(x-mu)^2)/(2*sigma^2))
(读作“给定 x、sigma 和 mu 的正态分布 n 由...给出”)到目前为止，很好；这与您拥有的功能相匹配。请注意，这里唯一真正的变量是 x :其他三个参数对于任何特定的高斯都是固定的。

现在来看一个数学事实:可以证明，所有高斯曲线都具有相同的形状，只是稍微移动了一点。所以我们可以使用 N(x|0,1,1) ，称为“标准正态分布”，只是将我们的结果转换回一般的高斯曲线。所以如果你有 N(x|0,1,1) 的积分，您可以轻松计算任何高斯的积分。这个积分出现得如此频繁以至于它有一个特殊的名字:误差函数erf .由于一些旧约定，它不完全是 erf ;还有一些加法和乘法因素也被随身携带。

如 Phi(z) = integral(N(x|0,1,1), -inf, z) ;也就是说，Phi(z)是从负无穷大到 z 的标准正态分布的积分，那么根据误差函数的定义是正确的
Phi(z) = 0.5 + 0.5 * erf(z / sqrt(2)) .

同样，如果 Phi(z | mu, sigma, n) = integral( N(x|sigma, mu, n), -inf, z) ;也就是说，Phi(z | mu, sigma, n)是给定参数的正态分布的积分 mu , sigma , 和 n从负无穷大到 z ，那么根据误差函数的定义是正确的
Phi(z | mu, sigma, n) = (n/2) * (1 + erf((x - mu) / (sigma * sqrt(2)))) .

看看the Wikipedia article on the normal CDF如果你想要更多细节或这个事实的证明。

好的，这应该是足够的背景解释。回到你的(编辑过的)帖子。你说“scipy.special 中的 erf(z) 需要我准确定义 t 最初是什么”。我不知道你的意思是什么；在哪里t (时间？)进入这个？希望上面的解释使错误函数变得神秘一点，现在更清楚为什么错误函数是工作的正确函数。

你的 Python 代码没问题，但我更喜欢闭包而不是 lambda:

def make_gauss(N, sigma, mu):
    k = N / (sigma * math.sqrt(2*math.pi))
    s = -1.0 / (2 * sigma * sigma)
    def f(x):
        return k * math.exp(s * (x - mu)*(x - mu))
    return f

使用闭包可以预先计算常量 k和 s ，因此每次调用返回的函数时都需要做更少的工作(如果您正在集成它，这可能很重要，这意味着它将被多次调用)。此外，我避免使用任何取幂运算符 ** ，这比写出平方要慢，并且将除法从内循环中取出并用乘法代替。我完全没有研究过它们在 Python 中的实现，但是从我上次使用原始 x87 程序集调整内部循环以获得纯速度时，我似乎记得加法、减法或乘法每个需要大约 4 个 CPU 周期，除以大约36 和大约 200 的幂。那是几年前的事了，所以对这些数字持保留态度；尽管如此，它还是说明了它们的相对复杂性。同样，计算 exp(x)蛮力方式是一个非常糟糕的主意；在编写 exp(x) 的良好实现时，您可以采取一些技巧。这使它比一般的 a**b 更快、更准确样式求幂。

我从未使用过常量 pi 和 e 的 numpy 版本；我一直坚持使用普通的旧数学模块版本。我不知道为什么你可能更喜欢任何一个。

我不确定您要使用 quad() 做什么称呼。 quad(gen_gauss, -inf, inf, (10,2,0))应该将重归一化的高斯从负无穷大积分到正无穷大，并且应该始终吐出 10(您的归一化因子)，因为高斯在实线上积分为 1。任何远离 10 的答案(我不希望正好是 10，因为 quad() 毕竟只是一个近似值)意味着某些地方出了问题......工作原理 quad() .

希望这已经揭开了一些困惑的神秘面纱，并解释了为什么错误函数是您问题的正确答案，以及如果您很好奇，如何自己做这一切。如果我的任何解释不清楚，我建议先快速浏览一下维基百科；如果您仍有疑问，请随时提出。