Modo migliore per scrivere una funzione Python che integra una gaussiana?

Question

Nel tentativo di utilizzare scipy quad metodo per integrare una gaussiana (diciamo che c'è una gaussiana metodo di gauss), ho avuto problemi passando i parametri necessari a gauss e lasciando quad per fare l'integrazione oltre la variabile corretta.Qualcuno ha un buon esempio di come utilizzare quad w/ multidimensionale funzione?

Ma questo mi ha portato a una più grande domanda circa il modo migliore per integrare una gaussiana in generale.Non riesco a trovare un gaussiana integrare in scipy (a mia sorpresa).Il mio piano era quello di scrivere una semplice funzione gaussiana e passare a quad (o forse ora una larghezza fissa integrator).Cosa vorresti fare?

Edit:A larghezza fissa che significa qualcosa come trapz che utilizza un fisso dx per calcolare le aree sotto la curva.

Che cosa sono venuto finora è un metodo di fare___gauss che restituisce una funzione lambda che possono andare in quad.In questo modo riesco a fare una funzione normale con media e varianza di cui ho bisogno prima di integrare.

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 ho provato il passaggio di una generale funzione gaussiana (che deve essere chiamato con x, N, mu, sigma) e la compilazione di alcuni dei valori mediante quad come

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

i parametri 10, 2 e 0 NON necessariamente corrispondono a N=10, sigma=2, v=0, che ha indotto la più estesa definizione.

Erf(z) in scipy.di speciale mi richiede di definire esattamente che cosa t viene inizialmente, ma è bello sapere che c'è.

Answer 1

Va bene, lei sembra essere piuttosto confuso su diverse cose.Cominciamo dall'inizio:lei ha parlato di un "multidimensionale funzione", ma poi andare a discutere di solito una variabile Gaussiana curva.Questo è non multidimensionale funzione:quando si integrano, è l'integrazione di una variabile (x).La distinzione è importante, perché ci è un mostro chiamato "distribuzione Gaussiana multivariata", che è un vero multidimensionale funzione e, se integrato, richiede l'integrazione di due o più variabili (che utilizza il costoso tecnica Monte Carlo ho accennato prima).Ma ti sembra di essere solo parlando di regolare una variabile Gaussiana, che è molto più facile lavorare con, integrare, e tutto ciò che.

Una variabile Gaussiana di distribuzione ha due parametri, sigma e mu, ed è una funzione di una sola variabile, ti indicano x.È, inoltre, sembrano essere portare in giro un parametro di normalizzazione n (il che è utile in un paio di applicazioni).La normalizzazione dei parametri sono di solito non incluso nei calcoli, dal momento che si può solo virare nuovamente alla fine (ricordate, l'integrazione è un operatore lineare: int(n*f(x), x) = n*int(f(x), x) ).Ma siamo in grado di portare in giro, se ti piace;la notazione mi piace per una distribuzione normale è quindi

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

(leggi come "la distribuzione normale di x dato sigma, mu, e n è dato da...") Così lontano, così buono;questo corrisponde alla funzione che hai.Si noti che l'unica la vera variabile qui è x:i tre parametri sono fisso per qualsiasi particolare Gaussiana.

Ora, per una formula matematica:è dimostrabile vero che tutti Gaussiana curve hanno la stessa forma, sono solo spostato un po'.Così siamo in grado di lavorare con N(x|0,1,1), chiamato la "distribuzione normale", e a tradurre i nostri risultati generali curva Gaussiana.Quindi, se avete l'integrale di N(x|0,1,1), si può facilmente calcolare l'integrale di qualsiasi Gaussiana.Questo integrale appare così spesso che ha un nome speciale:il la funzione di errore erf.A causa di alcune vecchie convenzioni, non esattamente erf;ci sono un paio di additivo e moltiplicativo fattori anche essere portato in giro.

Se Phi(z) = integral(N(x|0,1,1), -inf, z);che è, Phi(z) è l'integrale della distribuzione normale standard da meno infinito fino a z, poi è vero per la definizione della funzione di errore che

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

Allo stesso modo, se Phi(z | mu, sigma, n) = integral( N(x|sigma, mu, n), -inf, z);che è, Phi(z | mu, sigma, n) è l'integrale della distribuzione normale dato parametri mu, sigma, e n da meno infinito fino a z, poi è vero per la definizione della funzione di errore che

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

Guarda l'articolo di Wikipedia sul CDF normale se si desidera ulteriori informazioni o una prova di questo fatto.

Ok, questo dovrebbe essere sufficiente sfondo, spiegazione.Torna al tuo (a cura di) post.Si dice "erf(z) in scipy.di speciale mi richiede di definire esattamente che cosa t è inizialmente".Non ho idea di cosa si intende con questo;da dove viene t (tempo?) entra in tutto questo?Speriamo che la spiegazione di cui sopra ha demistificato la funzione di errore un po ' ed è più chiaro ora perché la funzione di errore è la funzione giusta per il lavoro.

Il codice Python è OK, ma io preferirei una chiusura su un 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

Utilizzando un chiusura consente precomputation di costanti k e s, così restituita la funzione dovrà fare meno lavoro ogni volta che si chiama (che può essere importante se si sta integrando, il che significa che sarà chiamato molte volte).Inoltre, ho evitato qualsiasi uso dell'operatore di elevamento a potenza **, che è più lento di scrivere la quadratura fuori, e issato la divisione uscita del loop interno e lo ha sostituito con un moltiplicarsi.Non ho guardato tutte le loro implementazione in Python, ma dalla mia ultima volta messa a punto di un ciclo interno per la velocità pura, utilizzo di materie prime 87 x il montaggio, mi sembra di ricordare che aggiunge, sottrae, o moltiplica circa 4 cicli di CPU per ogni, si divide su 36, e l'elevamento a potenza di circa 200.Che è stato un paio di anni fa, in modo da prendere quei numeri con un grano di sale;ancora, illustra la loro relativa complessità.Così, il calcolo exp(x) la forza bruta non è una cattiva idea,ci sono trucchi che si possono adottare durante la scrittura di una buona implementazione di exp(x) che rendono notevolmente più veloce e più preciso di un generale a**b stile di elevamento a potenza.

Non ho mai usato il numpy versione delle costanti pi ed e;Ho sempre bloccato con testo di matematica modulo versioni.Non so perché, si potrebbe preferire l'uno o l'altra.

Io non sono sicuro di quello che stai andando per il quad() chiamata. quad(gen_gauss, -inf, inf, (10,2,0)) dovrebbe integrare un rinormalizzati Gaussiana da meno infinito a più infinito, e deve sempre sputare 10 (il fattore di normalizzazione), dal momento che la Gaussiana si integra a 1 sulla retta reale.La risposta è lontano da 10 (non mi aspetto esattamente 10 dal quad() è solo un'approssimazione, dopo tutto), significa che qualcosa è sbagliato da qualche parte...difficile dire cosa c'è di sbagliato senza conoscere l'effettivo valore di ritorno e, eventualmente, il funzionamento interno di quad().

Speriamo che ha demistificato parte la confusione, e ha spiegato il perché della funzione di errore è la risposta giusta per il tuo problema, così come il modo di fare tutto da te, se siete curiosi.Se qualsiasi mia spiegazione non è stata chiara, io suggerisco di dare una rapida occhiata a Wikipedia prima;se avete ancora delle domande, non esitate a chiedere.

Answer 2

scipy navi con l'errore "funzione", aka Gaussiana integrale:

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

Answer 3

Presumo che tu stai gestione multivariata Gaussiane;se è così, SciPy ha già la funzione che si sta cercando:si chiama MVNDIST ("Distribuzione Normale Multivariata).Il SciPy documentazione è, come sempre, terribile, quindi non riesco nemmeno a trovare dove la funzione è sepolto, ma c'è da qualche parte.La documentazione è facilmente la parte peggiore di SciPy, e ha frustrato a me a non finire in passato.

Singole variabili Gaussiane basta usare il buon vecchio errore di funzione, di cui molte implementazioni sono disponibili.

Per affrontare il problema in generale, sì, come James Thompson cita, si vuole solo scrivere la tua distribuzione gaussiana funzione e dei mangimi per il quad().Se si può evitare la generalizzata integrazione, però, è una buona idea per fare in modo -- specializzato tecniche di integrazione per una particolare funzione (come MVNDIST usa) stanno per essere molto più veloce rispetto a una standard Monte Carlo multidimensionale di integrazione, che può essere estremamente lento per un'elevata precisione.

Answer 4

La distribuzione gaussiana è anche chiamato una distribuzione normale.Il cdf in scipy norma modulo fa quello che vuoi.

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

Perché non fare sempre il vostro integrazione da -infinito a +infinito, in modo da sapere sempre la risposta?(scherzo!)

La mia ipotesi è che l'unico motivo per cui non c'è già un fisso di funzione Gaussiana in SciPy è che si tratta di una banale funzione per scrivere.Il tuo suggerimento su come scrivere la vostra propria funzione e passando per quad integrare suoni eccellenti.Utilizza accettato SciPy strumento per fare questo, è il codice minimo sforzo per voi, ed è molto leggibile per altre persone anche se non hanno mai visto SciPy.

Cosa intendi esattamente per una larghezza fissa integratore?Vuoi dire che utilizza un algoritmo diverso rispetto a ciò che QUADPACK utilizza?

Edit:Per completezza, ecco qualcosa di simile a quello che mi piacerebbe provare per una Gaussiana con media 0 e deviazione standard 1, da 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 )

Che è un po 'brutta perché la funzione Gaussiana è un po' lungo, ma ancora abbastanza banale da scrivere.