编写一个Python函数，用于整合高斯分布的最佳方法是什么？

Question

编写一个Python函数，用于整合高斯分布的最佳方法是什么？

pythonscipygaussianintegral

7

在尝试使用scipy的quad方法来积分高斯函数（假设有一个名为gauss的高斯函数），我遇到了一些问题，无法将所需参数传递给gauss并使quad在正确的变量上进行积分。有人有使用带有多维函数的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是什么，但知道它在那里很好。

- vgm64

一个高斯分布的数字或数据。如果绘制出来，它看起来像一个隆起或“钟形曲线”。 - physicsmichael

1

口语化地讲，高斯通常用作一个名词来表示高斯曲线或分布（例如在维基百科条目中相当普遍）。我想我们也应该使用大写字母书写，但是SO是一个非常口语化的地方，不是吗？ - physicsmichael

5

“高斯分布是高斯分布，无论它修饰什么名词。请停止那些毫无意义的语义争论。” - temp2290

3

scipy.stats.norm.cdf 函数用于计算你的积分。 - John D. Cook

5个回答

17

scipy 自带了 "error function"，也称高斯积分：

import scipy.special
help(scipy.special.erf)

- Mr Fooz

3

您需要进行一次小的变量替换才能将误差函数（erf）转换为高斯累积分布函数（Gaussian CDF）。请参见此处的说明：http://www.johndcook.com/erf_and_normal_cdf.pdf。 - John D. Cook

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

- Chuck

3

为什么不总是从负无穷到正无穷进行积分，这样你就可以始终知道答案了呢？（开玩笑！）

我猜，SciPy中没有预设的高斯函数的唯一原因是它是一个微不足道的函数。您关于编写自己的函数并将其传递给QUAD以进行积分的建议听起来非常好。它使用了已接受的SciPy工具来执行此操作，对您来说，代码量很小，对其他人来说，即使他们从未见过SciPy，也很容易理解。

您所说的固定宽度积分器是什么意思？您是否意味着使用与QUADPACK使用不同的算法？

编辑：为了完整起见，以下是从0到+∞的均值为0，标准偏差为1的高斯函数的示例：

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 )

这有点丑，因为高斯函数有点长，但仍然很容易编写。

- James Thompson

3

我假设您正在处理多元高斯分布；如果是这样，那么SciPy已经有了您要找的函数：它被称为MVNDIST（“MultiVariate Normal DISTribution”）。 SciPy文档一如既往地糟糕，所以我甚至找不到该函数埋藏在哪里，但是它在里面。文档是SciPy最糟糕的部分，过去曾让我非常沮丧。

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

至于攻击问题的方法，正如James Thompson所提到的，您只需要编写自己的高斯分布函数并将其馈送给quad()。如果可以避免广义积分，则最好这样做--针对特定函数的专门积分技术（如MVNDIST使用的技术）比标准蒙特卡罗多维积分快得多，后者对于高精度可能非常慢。

- kquinn

网页内容由stack overflow 提供, 点击上面的

可以查看英文原文，
原文链接

- kquinn · Accepted Answer

好的，你似乎对几件事情感到困惑。让我们从头开始：你提到了“多维函数”，但接着又讨论了通常的单变量高斯曲线。这不是多维函数：当你对它进行积分时，你只积分一个变量（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))

(读作“给定 sigma、mu 和 n，x 的正态分布是...”) 到目前为止，一切顺利；这与您拥有的函数相匹配。请注意，这里唯一的真实变量是x：其他三个参数对于任何特定的高斯曲线都是固定的。

现在是一个数学事实：可以证明所有的高斯曲线具有相同的形状，只是稍微移动了一下。因此，我们可以使用N(x|0,1,1)，称为“标准正态分布”，并将结果转换回一般的高斯曲线。因此，如果你有N(x|0,1,1)的积分，你可以轻松地计算出任何高斯曲线的积分。这个积分出现得如此频繁，以至于它有一个特殊的名字：误差函数error functionerf。由于一些旧的惯例，它不是完全等于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))))。

如果您需要更多细节或证明这个事实，请参阅维基百科上的正态 CDF。

好的，这应该足够背景解释了。回到您（编辑后的）帖子。您说“scipy.special中的erf(z)需要我最初定义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 风格的指数运算更快、更准确。

我从未使用过 numpy 版本的常量 pi 和 e；我总是坚持使用普通的 math 模块版本。我不知道您可能更喜欢哪一个。

我不确定你使用quad()调用的目的是什么。quad(gen_gauss, -inf, inf, (10,2,0))应该从负无穷到正无穷积分一个经过重整化的高斯函数，应该总是输出10（您的归一化因子），因为高斯函数在实数线上积分为1。任何远离10的答案（毕竟quad()只是一个近似值）都意味着某些地方出了问题......很难说不知道实际返回值和可能的quad()内部工作情况。

希望这解开了一些困惑，并解释了为什么误差函数是您问题的正确答案，以及如果您好奇如何自行完成所有操作。如果我的解释有任何不清楚之处，建议先快速查看维基百科；如果您仍有疑问，请不要犹豫，随时提问。