6

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

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

编辑：固定宽度意味着类似于 trapz 的东西，它使用固定的 dx 来计算曲线下的面积。

到目前为止，我得到的是一个 make___gauss 方法，它返回一个 lambda 函数，然后可以进入四边形。这样我可以在积分之前用我需要的平均值和方差来制作一个正常的函数。

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(gen_gauss, -inf, inf, (10,2,0))

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

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

4

5 回答 5

34

好的，您似乎对几件事感到很困惑。让我们从头开始：您提到了一个“多维函数”，但接下来继续讨论通常的单变量高斯曲线。这不是一个多维函数：当你对它进行积分时，你只积分了一个变量 (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)的正态分布的积分，那么根据误差函数的定义，它是真的musigmanz

Phi(z | mu, sigma, n) = (n/2) * (1 + erf((x - mu) / (sigma * sqrt(2)))).

如果您想了解更多细节或证明这一事实，请查看有关普通 CDF 的 Wikipedia 文章。

好的，这应该是足够的背景解释。回到你的（编辑的）帖子。您说“ 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()。

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

于 2009-02-04T12:32:51.073 回答
15

scipy 附带“误差函数”，即高斯积分：

import scipy.special
help(scipy.special.erf)
于 2009-02-04T04:11:03.257 回答
4

高斯分布也称为正态分布。scipy norm 模块中的 cdf 函数可以满足您的需求。

from scipy.stats import norm
print norm.cdf(0.0)
>>>0.5

http://docs.scipy.org/doc/scipy/reference/generated/scipy.stats.norm.html#scipy.stats.norm

于 2011-01-18T13:49:59.387 回答
3

我假设您正在处理多元高斯；如果是这样，SciPy 已经拥有您正在寻找的功能：它被称为 MVNDIST（“MultiVariate Normal DISTribution”）。SciPy 文档一如既往地糟糕，所以我什至无法找到该函数的埋藏位置，但它在在某处。文档很容易成为 SciPy 中最糟糕的部分，过去一直让我感到沮丧。

单变量高斯仅使用良好的旧误差函数，其中有许多实现可用。

至于解决一般问题，是的，正如 James Thompson 所提到的，您只想编写自己的高斯分布函数并将其提供给 quad()。但是，如果您可以避免广义积分，那么这样做是个好主意——针对特定函数（如 MVNDIST 使用）的专用积分技术将比标准的蒙特卡罗多维积分快得多，后者可能非常慢为高精度。

于 2009-02-04T04:10:36.423 回答
2

为什么不总是进行从 -infinity 到 +infinity 的整合，以便您始终知道答案？（开玩笑！）

我的猜测是，SciPy 中还没有固定高斯函数的唯一原因是它是一个简单的函数。您关于编写自己的函数并将其传递给 quad 以进行集成的建议听起来很棒。它使用公认的 SciPy 工具来执行此操作，它对您来说是最少的代码工作，并且即使其他人从未见过 SciPy，它也非常易读。

固定宽度的积分器到底是什么意思？您的意思是使用与 QUADPACK 使用的算法不同的算法吗？

编辑：为了完整起见，这类似于我尝试使用平均值为 0 且标准差为 1 从 0 到 + 无穷大的高斯：

from scipy.integrate import quad
from math import pi, exp
mean = 0
sd   = 1
quad(lambda x: 1 / ( sd * ( 2 * pi ) ** 0.5 ) * exp( x ** 2 / (-2 * sd ** 2) ), 0, inf )

这有点难看，因为高斯函数有点长，但写起来还是很简单的。

于 2009-02-04T03:59:32.280 回答