首页 文章

编写集成高斯函数的Python函数的最佳方法是什么？

提问于
浏览 次
6

在尝试使用scipy的四边形方法来集成高斯（假设有一种名为gauss的高斯方法）时，我遇到了将所需参数传递给高斯并留下四边形来对正确变量进行积分的问题 . 有没有人有一个如何使用四维w /多维函数的好例子？

但这让我想到了一个关于整合高斯的最佳方法的更大问题 . 我没有在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(gen_gauss, -inf, inf, (10,2,0))

参数10,2和0不一定匹配N = 10，sigma = 2，mu = 0，这促使更多的扩展定义 .

scipy.special中的erf（z）将要求我确切地定义最初的t，但很高兴知道它就在那里 .

python scipy gaussian integral

5 回答

  • 12

    scipy船只带有“误差函数”，又称高斯积分：

    import scipy.special
help(scipy.special.erf)
    回复于
  • 3

    我假设你正在寻找：它甚至可以找到埋藏功能的地方，但是it's in there somewhere . 文档很容易成为SciPy中最糟糕的部分，并且让我在过去一直没有尽头 .

    单变量高斯只使用好的旧错误函数，其中有许多实现可用 .

    至于攻击问题一般，是的，正如James Thompson所提到的，你只想编写自己的高斯分布函数并将其提供给quad（） . 但是，如果可以避免广义集成，那么这样做是个好主意 - 特定函数的专用集成技术（如MVNDIST使用）将比标准的蒙特卡罗多维集成快得多，这可能非常慢高精度 .

    回复于
  • 3

    高斯分布也称为正态分布 . 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

    回复于
  • 32

    好吧，你似乎对几件事感到很困惑 . 让我们从头开始：你提到了"multidimensional function"，然后继续讨论通常的单变量高斯曲线 . 这不是一个多维函数：当你集成它时，你只集成一个变量（x） . 区别是很重要的，因为有一个名为"multivariate Gaussian distribution"的怪物，它是一个真正的多维函数，如果是集成的，需要整合两个或多个变量（使用我之前提到的昂贵的蒙特卡罗技术） . 但你似乎只是在谈论常规的单变量Gaussian，它更容易使用，集成，以及所有这些 .

    单变量高斯分布有两个参数 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 的正态分布 sigma ， mu ， n 由...给出”）到目前为止，这么好;这符合你所拥有的功能 . 请注意，这里唯一的真实变量是 x ：其他三个参数对于任何特定的高斯都是固定的 .

    现在有一个数学事实：可以证明所有高斯曲线具有相同的形状，它们只是稍微移动了一点 . 所以我们可以使用 N(x|0,1,1) ，称为"standard normal distribution"，并将结果转换回一般的高斯曲线 . 因此，如果你有 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 .

    好的，这应该是足够的背景解释 . 回到你的（编辑过的）帖子 . 你说"The erf(z) in scipy.special would require me to define exactly what t is initially" . 我不明白你的意思; 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) 的良好实现时，您可以采取一些技巧，使其比一般的 a**b 样式取幂更快，更准确 .

    我从未使用常数pi和e的numpy版本;我一直坚持使用普通的旧数学模块的版本 . 我不知道你为什么喜欢这两个 .

    我'm not sure what you'正在接受 quad() 电话 . quad(gen_gauss, -inf, inf, (10,2,0)) 应该将重正化的高斯从负无穷大加到无穷大，并且应该总是吐出10（你的归一化因子），因为高斯在实线上积分为1 . 任何远离10的答案（我不会指望自_1841809之后的10个只是近似值，毕竟）意味着某些东西被搞砸了......很难说在没有知道实际回报值和可能的内部运作的情况下搞砸了什么 quad() .

    希望这已经揭开了一些混乱的神秘面纱，并解释了为什么错误功能是您问题的正确答案，以及如果您感到好奇，如何自己完成 . 如果我的解释不清楚，我建议先快速浏览一下维基百科;如果您还有疑问，请不要犹豫 .

    回复于
  • 2

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

    我的猜测是，在SciPy中还没有高斯函数的唯一原因是它是一个很简单的函数 . 关于编写自己的函数并将其传递给四元组以集成声音的建议非常好 . 它使用可接受的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 )

    这有点难看，因为高斯函数有点长，但写起来仍然很简单 .

    回复于

相关问题