Melhor maneira de escrever uma função Python que integra um gaussian?

Question

Na tentativa de método quad uso de scipy para integrar um gaussian (permite dizer que não há um método de Gauss chamado gauss), eu estava tendo problemas para passar os parâmetros necessários para Gauss e deixando quad para fazer a integração sobre a variável correta. Alguém tem um bom exemplo de como usar quad w / a função multidimensional?

Mas isso me levou a um mais grande pergunta sobre a melhor maneira de integrar um gaussian em geral. I não encontrou um gaussian integrar em scipy (para minha surpresa). Meu plano era escrever uma função gaussiana simples e passá-lo para quad (ou talvez agora um integrador de largura fixa). O que você faria?

Edit:. Fixed-largura significando algo como Trapz que usa um dx fixo para áreas calcular sob uma curva

O que eu vim até agora é um método make___gauss que retorna uma função lambda que pode, então, entrar em quad. Dessa forma eu posso fazer uma função normal com a necessidade média e variância I 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)

Quando tentou passando uma função gaussiana geral (que tem de ser chamado com X, n, mu, e sigma) e o enchimento em alguns dos valores usando como quad

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

os parâmetros 10, 2, 0 e não coincidem necessariamente com N = 10, sigma = 2, mu = 0, o que levou a definição mais prolongado.

O erf (z) em scipy.special exigiria me definir exatamente o que t é inicialmente, mas bom saber que ele está lá.

Answer 1

Ok, você parece estar bastante confuso sobre várias coisas. Vamos começar no início: você mencionou uma "função multidimensional", mas, em seguida, passar a discutir a uma variável curva de Gauss habitual. Esta é não uma função multidimensional: quando você integrá-lo, você só integrar uma variável (x). A distinção é importante para fazer, porque há é um monstro chamado de "distribuição de Gauss multivariada", que é uma verdadeira função multidimensional e, se integrado, requer a integração ao longo de duas ou mais variáveis ??(que usa o Monte caro técnica Carlo mencionei antes). Mas você parece apenas estar falando de regular Gaussian uma variável, que é muito mais fácil trabalhar com, integrar, e tudo isso.

A distribuição de Gauss de uma variável tem dois parâmetros, sigma e mu, e é uma função de uma única variável vamos denotar x. Você também parecem estar transportando cerca de um n parâmetro de normalização (que é útil em algumas aplicações). parâmetros de normalização são geralmente não incluídos nos cálculos, desde que você pode apenas alinhavar-los de volta no final (lembre-se, a integração é um operador linear: int(n*f(x), x) = n*int(f(x), x)). Mas podemos carregá-lo se você gosta; a notação I como para uma distribuição normal é então

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

(leia-se "a distribuição normal dos x dada sigma, mu e n é dada por ...") Até aqui, tudo bem; isso corresponde a função que você tem. Observe que a única true variável aqui é x:. Os outros três parâmetros são fixa para qualquer Gaussian especial

Agora, para um fato matemático: é comprovadamente verdade que todas as curvas de Gauss têm a mesma forma, eles estão apenas deslocado ao redor um pouco. Assim, podemos trabalhar com N(x|0,1,1), chamado de "distribuição normal padrão", e apenas traduzir nossos resultados de volta para a curva de Gauss geral. Então se você tem a integral de N(x|0,1,1), você pode trivialmente calcular a integral de qualquer Gaussian. Esta integral aparece com tanta freqüência que ele tem um nome especial: o função de erro erf. Por causa de algumas velhas convenções, não é exatamente erf; há um par aditivo e factores multiplicativos também ser realizada em torno de.

Se Phi(z) = integral(N(x|0,1,1), -inf, z); isto é, Phi(z) é a integral da distribuição normal padrão de menos infinito até z, então é verdade pela definição da função de erro que

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

Da mesma forma, se Phi(z | mu, sigma, n) = integral( N(x|sigma, mu, n), -inf, z); isto é, Phi(z | mu, sigma, n) é a integral da distribuição normal de determinados parâmetros mu, sigma e n de menos infinito até z, então é verdade pela definição da função de erro que

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

Dê uma olhada o artigo da Wikipedia sobre o CDF normal se você quiser mais detalhes ou uma prova deste fato.

Ok, isso deve ser suficiente fundo explicação. Voltar para o seu post (editado). Você diz "O erf (z) em scipy.special exigiria me definir exatamente o que t é inicialmente". Eu não tenho nenhuma idéia do que você quer dizer com isso; onde é que t (tempo?) entrar neste em tudo? Esperemos que a explicação acima tenha desmistificado a função de erro um pouco e é mais claro agora por que motivo a função de erro é a função certa para o trabalho.

O seu código Python é OK, mas eu preferiria um encerramento de mais de um 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 um fecho permite precomputation de constantes k e s, de modo que o wil função retornoul necessidade de fazer menos trabalho cada vez que ele é chamado (que pode ser importante se você estiver integrando-o, o que significa que vai ser chamado muitas vezes). Além disso, tenho evitado qualquer uso do ** operador de exponenciação, que é mais lento do que apenas escrever a quadratura, e içada a divisão fora do circuito interno e substituiu-o com uma multiplicação. Eu não olhei para todos em sua implementação em Python, mas a partir de minha última vez sintonia um loop interno para a velocidade pura usando x87 matéria-montagem, eu me lembro que soma, subtrai, ou multiplica levar cerca de 4 ciclos de CPU cada um, divide sobre 36, e exponenciação cerca de 200. Isso foi um par de anos atrás, por isso tome esses números com um grão de sal; Ainda assim, ele ilustra sua relativa complexidade. Como assim, calcular exp(x) a forma de força bruta é uma idéia muito ruim; existem truques que você pode tomar quando se escreve uma boa implementação de exp(x) que o tornam significativamente mais rápido e mais preciso do que a exponenciação geral estilo a**b.

Eu nunca usei a versão numpy do pi constantes e e; Eu sempre preso com versões a planície antiga do módulo de matemática. Eu não sei porque você pode preferir um ou outro.

Eu não sei o que você está indo para com a chamada quad(). quad(gen_gauss, -inf, inf, (10,2,0)) deve integrar um Gaussian renormalizada de menos infinito a mais infinito, e sempre deve cuspir 10 (o fator de normalização), já que os integra Gauss a 1 sobre a linha real. Qualquer resposta longe de 10 (eu não esperaria exatamente 10 desde quad() é apenas uma aproximação, afinal de contas) significa algo é asneira em algum lugar ... difícil dizer o que é asneira sem saber o retorno real valor e, possivelmente, o funcionamento interno do quad().

Esperamos que tenha desmistificado algumas das confusões, e explicou por que a função de erro é a resposta certa para o seu problema, bem como a forma de fazer tudo sozinho, se você está curioso. Se algum de minha explicação não estava claro, eu sugiro dar uma olhada rápida na Wikipedia em primeiro lugar; se você ainda tiver dúvidas, não hesite em perguntar.

Answer 2

navios SciPy com a "função de erro", integrante aka Gaussian:

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

Answer 3

Eu suponho que você está segurando Gaussianas multivariadas; em caso afirmativo, SciPy já tem a função que você está procurando: ele é chamado MVNDIST ( "normal multivariada Distribuição) documentação O SciPy é, como sempre, terrível, então eu não posso mesmo encontrar onde a função é enterrado, mas é lá em algum lugar . a documentação é facilmente o pior parte de SciPy, e tem me frustrado sem fim no passado.

Gaussians de uma variável basta usar o bom função de erro de idade, dos quais muitas implementações estão disponíveis.

Como para atacar o problema em geral, sim, como James Thompson menciona, você só quer escrever sua própria função distribuição de Gauss e alimentá-lo para quad (). Se você pode evitar a integração generalizada, porém, é uma boa idéia para fazê-lo - técnicas de integração especializada para uma função específica (como usos MVNDIST) vão ser muito mais rápido do que uma integração multidimensional padrão Monte Carlo, que pode ser extremamente lento de alta precisão.

Answer 4

A distribuição de Gauss é também chamado de uma distribuição normal. A função cdf no módulo norma scipy faz o que quiser.

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 que não apenas sempre fazer a sua integração de -infinity a + infinito, de modo que você sempre sabe a resposta? (Brincadeira!)

Meu palpite é que a única razão que já não é uma função gaussiana enlatados em SciPy é que é uma função trivial para escrever. Sua sugestão sobre como escrever sua própria função e passá-la para quad para integrar sons excelente. Ele usa a ferramenta SciPy aceito para fazer isso, é esforço mínimo de código para você, e é muito legível para outras pessoas, mesmo que eles nunca viram SciPy.

O que exatamente você quer dizer com um integrador de largura fixa? Quer dizer usando um algoritmo diferente do que o que quer que QUADPACK está usando?

Edit: Para completar, algo está aqui como o que eu ia tentar por um Gaussian com a média 0 e desvio padrão 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 )

Isso é um pouco feio, porque a função de Gauss é um pouco longo, mas ainda bastante trivial para escrever.