интеграция функции 2d gaussian (python)

Question

Я определил 2D-гауссовский (без корреляции между независимыми переменными) с использованием параметров Area, sigmax и sigmay. Когда я выполняю интеграцию с (-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 

Выход

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 

Любая идея, почему это происходит? Я предполагаю, что это должно быть что-то простое, но после долгих часов смотреть на него я не вижу ничего плохого.

Answer 1

Правильный способ интеграции гауссовой ядра является использование Гаусса-Эрмита квадратур см here

Он реализован в Python как модуль SciPy http://docs.scipy.org/doc/numpy-1.10.1/reference/generated/numpy.polynomial.hermite.hermgauss.html

код, интеграл по гауссовой ядра должно быть & # 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 случай

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

Я думаю, что ваш гауссовский с 0,002 сигма слишком высок для квадратуры: Scipy игнорирует этот очень маленький пик и видит только нули повсюду. У Вас есть 2 решения:

• перенормировать функцию: ∫ ^б е ( х) ах = σ ∫ _а/σ^б/σ f ( u) du

• разрезать интеграл на многие части. Ниже приведен пример, который вычисляет интегралы от -бесконечности до -4 * сигма, то от -4 * сигма до 4 * сигма, а затем из 4 * сигма до бесконечности:

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 

я получить этот вывод :

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