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

Question

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

Pero esto me llevó a una pregunta más amplia 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:Ancho fijo significa algo así como trapz que usa un dx fijo para calcular áreas bajo una curva.

A lo que he llegado hasta ahora es a un método make___gauss que devuelve una función lambda que luego puede convertirse 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 debe llamarse con x, N, mu y sigma) y completar algunos de los valores usando quad como

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

los parámetros 10, 2 y 0 NO coincidían necesariamente con N=10, sigma=2, mu=0, lo que provocó una definición más amplia.

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

Answer 1

Bien, pareces estar bastante confundido acerca de varias cosas.Empecemos desde el principio:Mencionaste una "función multidimensional", pero luego pasas a discutir la curva gaussiana habitual de una variable.Esto es no una función multidimensional:cuando lo integras, solo integras una variable (x).Es importante hacer la distinción, porque hay es un monstruo llamado "distribución gaussiana multivariada" que es una verdadera función multidimensional y, si se integra, requiere la integración de dos o más variables (que utiliza la costosa técnica de Monte Carlo que mencioné antes).Pero parece que solo estás hablando del gaussiano regular de una variable, con el que es mucho más fácil trabajar, integrar y todo eso.

La distribución gaussiana de una variable tiene dos parámetros, sigma y mu, y es una función de una única variable que denotaremos x.También parece que llevas contigo un parámetro de normalización. n (que es útil en un par de aplicaciones).Los parámetros de normalización suelen ser no incluidos en los cálculos, ya que puedes volver a agregarlos al final (recuerda, la integración es un operador lineal: int(n*f(x), x) = n*int(f(x), x) ).Pero podemos llevarlo consigo si lo desea;la notación que 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))

(léase esto como "la distribución normal de x dado sigma, mu, y n está dado por...") Hasta ahora, todo bien;esto coincide con la función que tienes.Observe que el único verdadera variable aquí está x:los otros tres parámetros son fijado para cualquier gaussiano en particular.

Ahora un hecho matemático:Es demostrablemente cierto que todas las curvas gaussianas tienen la misma forma, sólo que están desplazadas un poco.Entonces podemos trabajar con N(x|0,1,1), llamada "distribución normal estándar", y simplemente traducimos nuestros resultados a la curva gaussiana general.Entonces si tienes la integral de N(x|0,1,1), puedes calcular trivialmente la integral de cualquier gaussiano.Esta integral aparece con tanta frecuencia que tiene un nombre especial:el función de error erf.Debido a algunas viejas convenciones, no es exactamente erf;También se llevan a cabo un par de factores aditivos y multiplicativos.

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

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

Asimismo, si Phi(z | mu, sigma, n) = integral( N(x|sigma, mu, n), -inf, z);eso es, Phi(z | mu, sigma, n) es la integral de la distribución normal dados los parámetros mu, sigma, y n desde menos infinito hasta z, entonces es cierto según la definición de la función de error que

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

Echa un vistazo a el artículo de Wikipedia sobre el CDF normal si quieres más detalles o una prueba de este hecho.

Bien, esa debería ser suficiente explicación de fondo.Regrese a su publicación (editada).Dices "El erf(z) en scipy.special requeriría que defina exactamente qué es t inicialmente".No tengo idea de lo que quieres decir con esto;donde hace t (¿tiempo?) entrar en esto?Esperemos que la explicación anterior haya desmitificado un poco la función de error y ahora quede 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 preferiría un cierre a 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

El uso de un cierre permite el cálculo previo de constantes k y s, por lo que la función devuelta necesitará hacer menos trabajo cada vez que se llame (lo cual puede ser importante si la estás integrando, lo que significa que se llamará muchas veces).Además, he evitado cualquier uso del operador de exponenciación. **, que es más lento que simplemente escribir la cuadratura, y sacó la división del bucle interno y la reemplazó con una multiplicación.No he examinado en absoluto su implementación en Python, pero desde la última vez que ajusté un bucle interno para obtener velocidad pura usando un ensamblaje x87 sin formato, creo recordar que las sumas, restas o multiplicaciones requieren aproximadamente 4 ciclos de CPU cada una, y divide aproximadamente 36 y exponenciación alrededor de 200.Eso fue hace un par de años, así que tomen esas cifras con cautela;aún así, ilustra su relativa complejidad.Así mismo, calculando exp(x) el método de la fuerza bruta es una muy mala idea;Hay trucos que puedes seguir al escribir una buena implementación de exp(x) que lo hacen significativamente más rápido y más preciso que un general a**b exponenciación del estilo.

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

No estoy seguro de qué estás buscando con el quad() llamar. quad(gen_gauss, -inf, inf, (10,2,0)) Debería integrar un gaussiano renormalizado desde menos infinito hasta más infinito, y siempre debería escupir 10 (su factor de normalización), ya que el gaussiano se integra a 1 sobre la línea real.Cualquier respuesta lejos de 10 (no esperaría exactamente 10 desde quad() es sólo una aproximación, después de todo) significa que algo está mal en alguna parte...Es difícil decir qué es lo que está mal sin conocer el valor de retorno real y posiblemente el funcionamiento interno de quad().

Con suerte, esto ha desmitificado parte de la confusión y ha explicado por qué la función de error es la respuesta correcta a su problema, además de cómo hacerlo todo usted mismo si tiene curiosidad.Si alguna de mis explicaciones no quedó clara, sugiero que primero eche un vistazo rápido a Wikipedia;Si todavía tienes preguntas, no dudes en preguntar.

Answer 2

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

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

Answer 3

Asumo que está manejando gaussianas multivariantes; si es así, SciPy ya cuenta con la función que está buscando: se llama MVNDIST ( "distribución normal multivariante) La documentación SciPy es, como siempre, terrible, así que ni siquiera puede encontrar dónde está enterrada la función, pero que está en alguna parte . la documentación es fácilmente la peor parte de SciPy, y me ha frustrado a ningún extremo en el pasado.

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

En cuanto a atacar el problema en general, sí, como se menciona a James Thompson, lo que desea escribir su propia función de distribución gaussiana y alimentar a quad (). Si se puede evitar la integración generalizada, sin embargo, es una buena idea hacerlo - técnicas de integración especializados para una función particular (como utiliza MVNDIST) van a ser mucho más rápido que una integración multidimensional estándar de Monte Carlo, que puede ser extremadamente lenta para alta precisión.

Answer 4

La distribución gaussiana también se llama una distribución normal. La función cdf en el módulo 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

Answer 5

¿Por qué no hacer siempre su integración desde -infinity a + infinito, de modo que siempre sabe la respuesta? (Broma!)

Mi conjetura es que la única razón de que no hay ya una función gaussiana enlatado en SciPy es que es una función trivial para escribir. Su sugerencia sobre cómo escribir su propia función y pasarla a quad integrar suena excelente. Utiliza la herramienta SciPy aceptado para hacer esto, es posible código mínimo para usted, y es muy fácil de leer para otras personas, incluso si nunca han visto SciPy.

¿Qué es exactamente lo que quiere decir con un integrador de ancho fijo? Qué quiere decir el uso de un algoritmo diferente a lo QUADPACK está utilizando?

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

Eso es un poco feo porque la función de Gauss es un poco largo, pero sigue siendo bastante trivial para escribir.