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

在尝试使用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,但很高兴知道它就在那里.

1> kquinn..：

---

好吧,你似乎对几件事感到很困惑.让我们从头开始:你提到了一个"多维函数",然后继续讨论通常的单变量高斯曲线.这不是一个多维函数:当你集成它时,你只集成一个变量(x).这种区别是重要的是要的,因为有是一个称为"多元高斯分布",这是一个真正的多维函数的怪物,如果集成,需要在两个或多个变量集成(使用昂贵的蒙特卡洛技术我前面提到过).但你似乎只是在谈论常规的单变量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))

(读取为"正态分布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)))).

如果您想了解更多细节或证明这一事实,请查看普通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().

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

---

@ vgm64:实际上,erf可以正常工作:假设您想要从mu-delta集成到mu + delta.积分只是Phi(mu + delta | mu,sigma,n) - Phi(mu-delta | mu,sigma,n):上面用erf()定义的Phi函数是高斯的反导数.

---

请修复`make_gauss`最后一行的缩进.

2> Mr Fooz..：

---

scipy船只带有"误差函数",又称高斯积分:

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

---

您需要对变量进行少量更改才能将erf转换为高斯CDF.请参阅此处的说明:http://www.johndcook.com/erf_and_normal_cdf.pdf