二维高斯函数（python）的积分

Question

我用面积、sigmax和sigmay参数定义了2D高斯（独立变量之间没有相关性）。 当我在两个变量中从（-inf，inf）进行积分时，我只得到sigmax和sigmay为1时的面积。在

import numpy as np
import scipy.integrate as sci

class BeamDistribution(object):
    def __init__(self, Ipeak, sigmax, sigmay):
        print Ipeak, sigmax, sigmay
        self.__Ipeak = Ipeak
        self.__sigmax = sigmax
        self.__sigmay = sigmay

    def value(self, x, y):
        factor = self.__Ipeak/(2.*np.pi*self.__sigmax * self.__sigmay)
        factorx = np.exp(-x**2/(2.*self.__sigmax**2))
        factory = np.exp(-y**2/(2.*self.__sigmay**2))
        return factor*factorx*factory

    def integral(self, a, b, c, d):
        integration = sci.dblquad(self.value, a, b, lambda x: c, lambda x: d,
                                  epsrel = 1e-9, epsabs = 0)
#        sci.quad_explain()
        return integration

    def __call__(self, x, y):
        return self.value(x, y)


if __name__ == "__main__":
    Ipeak = 65.0e-3
    sigmax = 0.2e-3
    sigmay = 0.3e-3
    limit = np.inf
    my_beam_class = BeamDistribution(Ipeak, sigmax, sigmay)
    total = my_beam_class.integral(-limit, limit, -limit, limit)
    print "Integrated total current ",total," of Ipeak ", Ipeak

    my_beam_class = BeamDistribution(Ipeak, 1, 1)
    total = my_beam_class.integral(-limit, limit, -limit, limit)
    print "Integrated total current ",total," of Ipeak ", Ipeak

输出是

^{pr2}$

你知道为什么会这样吗？我想这应该是简单的事情，但经过几个小时的观察，我看不出有什么不对的。在

Answer 1

积分高斯核的正确方法是使用高斯-厄米积分，见here

它在Python中实现为SciPy模块http://docs.scipy.org/doc/numpy-1.10.1/reference/generated/numpy.polynomial.hermite.hermgauss.html

代码，高斯核上的积分应该是√（π）

import math
import numpy as np

a,w = np.polynomial.hermite.hermgauss(32)

print(a)
print(w)

def f(x):
    return 1.0

s = 0.0
for k in range(0,len(a)):
    s += w[k]*f(a[k])

print(s - math.sqrt(math.pi))

二维案例

^{pr2}$

Answer 2

我认为0.002西格玛的高斯函数对于求积来说太过尖峰了：Scipy忽略了这个非常小的峰值，到处都是零。您有两种解决方案：

• σ函数

• 把积分分成许多部分。下面是一个计算从-无穷到-4*sigma，然后从-4*sigma到4*sigma，再从4*sigma到无穷大的积分：

def integral(self, a, b, c, d):
    integration =0
    nsigmas=4
    for intervalx in [(a,-nsigmas*sigmax),(-nsigmas*sigmax,nsigmas*sigmax),(nsigmas*sigmax,b)]:
        for intervaly in [(c,-nsigmas*sigmay),(-nsigmas*sigmay,nsigmas*sigmay),(nsigmas*sigmay,d)]:
                integration+= sci.dblquad(self.value, intervalx[0], intervalx[1], lambda x: intervaly[0], lambda x: intervaly[1],
                              epsrel = 1e-9, epsabs = 0)[0]
    #        sci.quad_explain()
    return integration

我得到这个输出：

^{pr2}$