CODEWITHSUNDEEP
pythonscipygaussiannumerical-integration
zunbeltz
zunbeltz

Reputation: 289

integration of 2d gaussian function (python)

I have defined a 2D Gaussian (without correlation between the independent variables) using the Area, sigmax and sigmay parameters. When I do a integration from (-inf, inf) in both variables I only get the Area when sigmax and sigmay are 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

The output is

0.065 0.0002 0.0003
Integrated total current  (7.452488478001055e-32, 6.855160478762106e-41)  of Ipeak  0.065
0.065 1 1
Integrated total current  (0.4084070449667172, 1.0138233535120856e-11)  of Ipeak  0.065

Any idea why is this happening? I guess it should be something simple, but afters hours looking to it I can not see anything wrong.

Upvotes: 2

Views: 4750

Answers (2)

Severin Pappadeux
Severin Pappadeux

Reputation: 20130

Proper way to integrate gaussian kernels is to use Gauss-Hermite quadratures, see here

It is implemented in Python as SciPy module http://docs.scipy.org/doc/numpy-1.10.1/reference/generated/numpy.polynomial.hermite.hermgauss.html

Code, integral over gaussian kernel should be √(π) 

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))

2D case

Ipeak = 0.065
sigmax = 0.2e-3
sigmay = 0.3e-3
sqrt2 = math.sqrt(2.)

def h(x, y):
    return Ipeak*1.0/math.pi

s = 0.0
for k in range(0, len(a)):
    x = sqrt2 * sigmax * a[k]
    t = 0.0
    for l in range(0, len(a)):
        y = sqrt2 * sigmay * a[l]
        t += w[l] * h(x, y)
    s += w[k]*t

print(s)

Upvotes: 2

JPG
JPG

Reputation: 2264

I think your gaussian with 0.002 sigma is far too peaked for a quadrature: Scipy ignores this very little peak and sees only zeros everywhere. You have 2 solutions:

  • renormalize the function: ∫ab f(x) dx = σ ∫a/σb/σf(u) du

  • slice the integral in many pieces. Here is an example which computes integrals from -infinity to -4*sigma, then from -4*sigma to 4*sigma, and then from 4*sigma to infinity:

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

I obtain this output:

0.065 0.0002 0.0003
Integrated total current  0.06499999987174367  of Ipeak  0.065
0.065 1 1
Integrated total current  0.06500000000019715  of Ipeak  0.065

Upvotes: 1

Related Questions