¿La mejor manera de escribir una función de Python que integre un gaussiano?

Question

Al intentar usar el método cuádruple de scipy para integrar un gaussiano (digamos que hay un método gaussiano llamado gauss), tuve problemas para pasar los parámetros necesarios a gauss y dejar quad para hacer la integración sobre la variable correcta. ¿Alguien tiene un buen ejemplo de cómo usar quad w/a una función multidimensional?

Pero esto me llevó a una gran pregunta sobre la mejor manera de integrar un gaussiano en general. No encontré una integración gaussiana en scipy (para mi sorpresa). Mi plan era escribir una función gaussiana simple y pasarla a quad (o tal vez ahora a un integrador de ancho fijo). ¿Qué harías?

Editar: Anchura fija significa algo así como trapz que usa un dx fijo para calcular áreas bajo una curva.

Lo que he llegado hasta ahora es un método make___gauss que devuelve una función lambda que luego puede entrar en quad. De esta manera puedo hacer una función normal con el promedio y la varianza que necesito antes de integrar.

def make_gauss(N, sigma, mu): 
    return (lambda x: N/(sigma * (2*numpy.pi)**.5) * 
      numpy.e ** (-(x-mu)**2/(2 * sigma**2))) 

quad(make_gauss(N=10, sigma=2, mu=0), -inf, inf) 

Cuando intenté pasar una función gaussiana general (que necesita ser llamado con x, N, mu, y sigma) y el llenado en algunos de los valores utilizando quad como

quad(gen_gauss, -inf, inf, (10,2,0)) 

los parámetros 10, 2 y 0 NO necesariamente coinciden con N = 10, sigma = 2, mu = 0, lo que provocó la definición más extendida.

El erf (z) en scipy.special requeriría que defina exactamente qué t es inicialmente, pero es bueno saber que está allí.

Answer 1

Está bien, parecen ser bastante confundido acerca de varias cosas. Comencemos por el principio: mencionaste una "función multidimensional", pero luego pasas a discutir la curva gaussiana usual de una variable. Esto es no una función multidimensional: cuando la integra, solo integra una variable (x). Es importante hacer la distinción, porque es un monstruo llamado "distribución gaussiana multivariante" que es una verdadera función multidimensional y, si está integrada, requiere integrar más de dos o más variables (que usa la costosa técnica de Monte Carlo que mencioné antes)) Pero parece que solo estás hablando del Gaussiano regular de una variable, que es mucho más fácil de trabajar, integrar y todo eso.

La distribución de Gauss de una variable tiene dos parámetros, sigma y mu, y es una función de una sola variable que vamos a denotamos x. También parece llevar un parámetro de normalización n (que es útil en un par de aplicaciones). Los parámetros de normalización generalmente son no incluidos en los cálculos, ya que puede volverlos a conectar al final (recuerde, la integración es un operador lineal: int(n*f(x), x) = n*int(f(x), x)). Pero podemos llevarlo a cabo si lo desea; la notación me gusta para una distribución normal es entonces

N(x | mu, sigma, n) := (n/(sigma*sqrt(2*pi))) * exp((-(x-mu)^2)/(2*sigma^2))

(leído que como "la distribución normal de x dado sigma, mu y n está dada por ...") Hasta ahora, todo va bien; esto coincide con la función que tienes. Observe que la única verdadera variable aquí es x: los otros tres parámetros son fijo para cualquier gaussiano en particular.

Ahora, por un hecho matemático: es probable que todas las curvas de Gauss tengan la misma forma, simplemente se desplazan un poco. Entonces podemos trabajar con N(x|0,1,1), llamado "distribución normal estándar", y simplemente traducimos nuestros resultados a la curva general de Gauss. Entonces, si tiene la integral de N(x|0,1,1), puede calcular trivialmente la integral de cualquier gaussiano. Esta integral aparece con tanta frecuencia que tiene un nombre especial: la función de error erf. Debido a algunas convenciones antiguas, no es exactamenteerf; hay un par de factores aditivos y multiplicativos que también se transportan.

Si Phi(z) = integral(N(x|0,1,1), -inf, z); es decir, Phi(z) es la integral de la distribución normal estándar de menos infinito hasta z, entonces es cierto por la definición de la función de error que

Phi(z) = 0.5 + 0.5 * erf(z/sqrt(2)).

Del mismo modo, si Phi(z | mu, sigma, n) = integral(N(x|sigma, mu, n), -inf, z); es decir, Phi(z | mu, sigma, n) es la integral de los parámetros de distribución normal dadas mu, sigma, y n de menos infinito hasta z, entonces es cierto por la definición de la función de error que

Phi(z | mu, sigma, n) = (n/2) * (1 + erf((x - mu)/(sigma * sqrt(2)))).

Eche un vistazo a the Wikipedia article on the normal CDF si quiere más detalles o una prueba de esto.

Bueno, eso debería ser suficiente explicación de fondo. De vuelta a su publicación (editada). Usted dice "El erf (z) en scipy.special requeriría que defina exactamente qué t es inicialmente". No tengo idea de lo que quieres decir con esto; ¿dónde entra t (¿tiempo?) en esto? Afortunadamente, la explicación anterior ha desmitificado un poco la función de error y ahora está más claro por qué la función de error es la función correcta para el trabajo.

Su código Python está bien, pero yo preferiría un cierre sobre una 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 

Usando un cierre permite precomputation de constantes k y s, por lo que la función devuelta tendrá que hacer menos trabajo cada vez se llama (lo cual puede ser importante si lo está integrando, lo que significa que será llamado muchas veces). Además, he evitado el uso del operador de exponenciación **, que es más lento que solo escribir la cuadratura, y he levantado la división del ciclo interno y la he reemplazado por una multiplicación. No he visto en absoluto su implementación en Python, pero desde la última vez que sintonicé un bucle interno para velocidad pura usando el ensamblaje x87 sin procesar, me parece recordar que agrega, resta o multiplica toma alrededor de 4 ciclos de CPU cada uno, se divide sobre 36, y exponenciación de aproximadamente 200. Eso fue hace un par de años, así que tome esos números con un grano de sal; aún así, ilustra su complejidad relativa. Además, calcular exp(x) el modo de fuerza bruta es una muy mala idea; hay trucos que puede tomar al escribir una buena implementación de exp(x) que lo hacen significativamente más rápido y más preciso que una exponenciación de estilo general a**b.

Nunca he usado la versión numpy de las constantes pi y e; Siempre me he quedado con las versiones simples del módulo de matemáticas. No sé por qué podrías preferir cualquiera de los dos.

No estoy seguro de lo que está buscando con la llamada quad(). quad(gen_gauss, -inf, inf, (10,2,0)) debe integrar un gaussiano renormalizado de menos infinito a más infinito, y siempre debe escupir 10 (su factor de normalización), ya que Gaussian se integra a 1 sobre la línea real. Cualquier respuesta lejos de 10 (no esperaría que exactamente 10 ya que quad() es sólo una aproximación) significa que algo está arruinado en algún lado ... es difícil decir lo que se estropeó sin conocer el valor de retorno real y posiblemente el valor interno funcionamientos de quad().

Esperemos que esto haya desmitificado parte de la confusión y haya explicado por qué la función de error es la respuesta correcta a su problema, y ​​cómo hacerlo todo usted mismo si tiene curiosidad. Si alguna de mis explicaciones no estaba clara, sugiero echar un vistazo rápido a Wikipedia primero; si todavía tiene preguntas, no dude en preguntar.

Answer 2

¿Por qué no hacer siempre su integración de infinito a infinito, para que siempre sepa la respuesta? (¡en broma!)

Supongo que la única razón por la que ya no hay una función gaussiana enlatada en SciPy es que es una función trivial para escribir. Su sugerencia acerca de escribir su propia función y pasarla a quad para integrar sonidos excelentes. Utiliza la herramienta SciPy aceptada para hacer esto, es un esfuerzo de código mínimo para ti, y es muy legible para otras personas, incluso si nunca han visto SciPy.

¿Qué quiere decir exactamente con un integrador de ancho fijo? ¿Quiere decir usar un algoritmo diferente de lo que QUADPACK está usando?

Editar: Para completar, aquí hay algo parecido a lo que iba a tratar de una gaussiana con la media de 0 y una desviación estándar de 1 a 0 a + infinito:

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) 

que es un poco fea porque el Gauss la función es un poco larga, pero aún bastante trivial de escribir.

Answer 3

Supongo que está manejando Gaussianos multivariantes; si es así, SciPy ya tiene la función que estás buscando: se llama MVNDIST ("MultiVariate Normal DisTribution"). La documentación de SciPy es, como siempre, terrible, así que ni siquiera puedo encontrar dónde está enterrada la función, pero it's in there somewhere. la documentación es fácilmente la peor parte de SciPy, y me ha frustrado a ningún extremo en el pasado.

sola variable gaussianas sólo tiene que utilizar la buena función de error de edad, de las cuales muchas implementaciones están disponibles.

cuanto a atacando el problema en general, sí, como menciona James Thompson, solo quieres escribir tu propia función de distribución gaussiana y alimentarla a quad().Sin embargo, si puede evitar la integración generalizada, es una buena idea hacerlo: las técnicas de integración especializadas para una función en particular (como MVNDIST usa) van a ser mucho más rápidas que una integración multidimensional Monte Carlo estándar, que puede ser extremadamente lenta. para una alta precisión.

Answer 4

barcos SciPy con la "función de error", también conocido como Integral de Gauss:

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

Answer 5

La distribución gaussiana también se denomina distribución normal. La función cdf en el módulo de norma scipy hace lo que quiere.

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