Integration der 2d-Gauss-Funktion (Python)

Question

Ich habe einen 2D-Gaussian (ohne Korrelation zwischen den unabhängigen Variablen) mit den Parametern Area, Sigmax und Sigmay definiert. Wenn ich eine Integration von (inf, inf) in beiden Variablen ich die Gegend erhalten nur, wenn Sigmax und sigmay sind 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 

Der Ausgang ist

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 

Jede Idee, warum das so ist Ereignis? Ich denke, es sollte etwas Einfaches sein, aber nach Stunden, wenn ich darauf schaue, kann ich nichts falsches sehen.

Answer 1

Die richtige Art und Weise Gaußkerne zu integrieren ist Gauss-Hermite Quadraturen zu verwenden, finden Sie here

Es ist in Python als SciPy Modul implementiert ist http://docs.scipy.org/doc/numpy-1.10.1/reference/generated/numpy.polynomial.hermite.hermgauss.html

Kodex soll Integral über Gaußkerns sein & # x221a; (π)

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 Fall

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) 

Answer 2

Ich denke, dein Gaussian mit 0,002 Sigma ist viel zu hoch für eine Quadratur: Scipy ignoriert diesen sehr kleinen Peak und sieht überall nur Nullen. Sie haben 2 Lösungen:

• renormieren die Funktion: ∫ _a ^b f ( x) dx = σ ∫ _a/σ^b/σ f ( u) du

• schneiden Sie das Integral in vielen Stücken. Hier ist ein Beispiel, das von Integralen -Infinity berechnet auf -4 * Sigma, dann von -4 * sigma bis 4 * sigma, und dann von 4 * sigma bis unendlich:

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 diese Ausgabe erhalten, :

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