编写包含高斯函数的python函数的最佳方法?

好吧,你似乎对一些事情很困惑。让我们从一开始:您提到了一个“多维函数”,然后继续讨论通常的单变量高斯曲线。这是 不 多维函数:当你整合它时,你只整合一个变量(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 鉴于 西格玛 , 穆 和 n 到目前为止,很好,这与你所得到的函数相匹配。注意,只有 真变量 这里是 X :其他三个参数是 固定的 对于任何特定的高斯。

现在,对于一个数学事实:可以证明,所有的高斯曲线都有相同的形状,只是稍微移动了一点。所以我们可以合作 N(x|0,1,1) 称为“标准正态分布”,只需将我们的结果转换回一般的高斯曲线。所以如果你有 n(x0,1,1) ,你可以简单地计算任意高斯积分。这个积分出现得如此频繁,以至于有一个特殊的名称:the 误差函数 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) 是正态分布给定参数的积分 穆 , 西格玛 和 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中的实现,但是从我上次使用raw x87程序集调整内部循环以获得纯粹的速度开始,我似乎记得加法、减法或乘法每个大约需要4个CPU周期,除法约36,求幂约200。那是几年前的事了,所以用一粒盐来计算这些数字,它仍然说明了它们的相对复杂性。同时,计算 exp(x) 蛮力方法是一个非常糟糕的主意;在编写一个好的实现 EXP(X) 这使得它比一般的方法更快更准确。 a**b 样式求幂。

我从来没有使用过常量pi和e的麻木版本;我一直坚持使用简单的旧数学模块的版本。我不知道你为什么会喜欢这两个。

我不知道你要干什么 quad() 打电话。 quad(gen_gauss, -inf, inf, (10,2,0)) 应该将重正化高斯从负无穷大积分到正无穷大,并且应该总是吐出10(标准化因子),因为高斯在实线上积分到1。有10个以外的答案吗(我没想到 确切地 10以来 四() 毕竟,这只是一个近似值)意味着某些事情在某个地方被搞砸了…很难说什么是搞砸了,不知道实际的返回值,可能还有内部的工作 四() .

希望这已经解开了一些困惑,并解释了为什么错误函数是解决问题的正确答案,以及如果你好奇的话如何自己去做。如果我的解释不清楚,我建议你先快速浏览一下维基百科;如果你还有问题,不要犹豫。

scipy附带“误差函数”,即高斯积分:

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

我假设你在处理多变量高斯函数;如果是这样,那么scipy已经有了你想要的函数:它被称为mvndist(“多变量正态分布”)。scipy文档和以前一样糟糕,所以我甚至找不到函数的隐藏位置,但是 it's in there somewhere . 文档很容易是scipy最糟糕的部分,让我在过去一直感到沮丧。

单变量Gaussian只使用好的旧错误函数,其中有许多实现可用。

对于一般的攻击问题,是的,正如詹姆斯·汤普森所提到的,您只需要编写自己的高斯分布函数并将其馈送给Quad()。如果你能避免广义积分,那么这样做是一个好主意——特定函数的专门积分技术(如mvndist使用)将比标准的蒙特卡洛多维积分快得多,这对于高精度来说非常慢。

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

为什么不总是把你的积分从-无穷大到+无穷大,这样你就总能知道答案?(开玩笑!)

我的猜测是,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 )

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