La meilleure façon d'écrire une fonction Python qui intègre une gaussienne?

Question

En essayant d'utiliser la méthode de quad scipy pour intégrer une gaussienne (disons que il y a une méthode gaussienne gauss nom), j'avais des problèmes en passant les paramètres nécessaires pour gauss et laissant quad à faire l'intégration sur la variable correcte. Quelqu'un at-il un bon exemple de la façon d'utiliser quad w / une fonction multidimensionnelle?

Mais cela me conduit à une question plus grande sur la meilleure façon d'intégrer une gaussienne en général. Je n'ai pas trouvé une gaussienne à intégrer scipy (à ma grande surprise). Mon plan était d'écrire une simple fonction gaussienne et le transmettre à quad (ou peut-être maintenant un intégrateur de largeur fixe). Que feriez-vous?

Edit:. Sens de largeur fixe quelque chose comme trapz qui utilise un dx fixe pour calculer les zones sous la courbe

Ce que je suis venu à ce jour est une make___gauss méthode qui retourne une fonction lambda qui peut ensuite en quad. De cette façon, je peux faire une fonction normale avec la moyenne et la variance j'ai besoin avant d'intégrer.

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)

Quand j'ai essayé le passage d'une fonction gaussienne générale (qui doit être appelé avec x, N, mu et sigma) et remplissant certaines des valeurs en utilisant quad comme

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

les paramètres 10, 2 et 0 ne correspond pas nécessairement à N = 10, sigma = 2, mu = 0, ce qui a poussé la définition plus étendue.

Le erf (z) scipy.special me demander de définir exactement ce que t est au départ, mais agréable de savoir qu'il est là.

Answer 1

D'accord, vous semblez être assez confus au sujet de plusieurs choses. Commençons par le début: vous avez parlé d'une « fonction multidimensionnelle », mais aller à discuter de la courbe de Gauss à une variable habituelle. Ceci est pas une fonction multidimensionnelle: lorsque vous intégrez, vous n'intégrez une variable (x). La distinction est importante à faire, parce qu'il est un monstre appelé une « distribution gaussienne à plusieurs variables », qui est une véritable fonction multidimensionnelle et, si elle est intégrée, nécessite d'intégrer plus de deux ou plusieurs variables (qui utilise le cher Monte Carlo technique je l'ai mentionné auparavant). Mais vous semblez juste parler de la normale gaussienne à une variable, ce qui est beaucoup plus facile de travailler avec, d'intégrer et de tout cela.

La seule variable distribution gaussienne a deux paramètres, et sigma mu, et est fonction d'une seule variable nous noterons x. Vous semblez aussi être trimballer un paramètre de normalisation n (ce qui est utile dans quelques applications). paramètres sont généralement pas inclus dans les calculs, puisque vous pouvez simplement les virer de bord en arrière sur à la fin (rappelez-vous, l'intégration est un opérateur linéaire: int(n*f(x), x) = n*int(f(x), x)). Mais nous pouvons transporter si vous le souhaitez; la notation J'aime une distribution normale est alors

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

(lire que « la distribution normale des données N(x|0,1,1) erf, Phi(z) = integral(N(x|0,1,1), -inf, z) et est donnée par Phi(z) ... ») Jusqu'à présent, si bon; cela correspond à la fonction que vous avez. Notez que la seule true variable ici est z:. Les trois autres paramètres sont fixe pour tout gaussienne

particulier

Maintenant, pour un fait mathématique: il est prouvable vrai que toutes les courbes de Gauss ont la même forme, ils sont juste déplacés autour d'un peu. Ainsi, nous pouvons travailler avec Phi(z) = 0.5 + 0.5 * erf(z / sqrt(2)), appelé la « distribution normale », et juste traduire nos résultats à la courbe de Gauss générale. Donc, si vous avez l'intégrale de Phi(z | mu, sigma, n) = integral( N(x|sigma, mu, n), -inf, z), vous pouvez trivialement calculer l'intégrale de toute gaussienne. Cette intégrale apparaît si souvent qu'il a un nom spécial: fonction d'erreur Phi(z | mu, sigma, n). En raison de certaines vieilles conventions, ce n'est pas exactement Phi(z | mu, sigma, n) = (n/2) * (1 + erf((x - mu) / (sigma * sqrt(2)))); il y a un additif couple et les facteurs multiplicatifs également en cours autour.

Si t; à savoir, l'intégrale est k de la distribution normale de moins l'infini jusqu'à s, il est vrai par la définition de la fonction d'erreur

**.

De même, si exp(x); à savoir, l'intégrale est a**b de la distribution normale étant donné les paramètres quad(), quad(gen_gauss, -inf, inf, (10,2,0)) et de l'infini <=> moins jusqu'à <=>, il est vrai par la définition de la fonction d'erreur

<=>.

Jetez un oeil à l'article de Wikipédia sur la normale CDF si vous voulez plus de détails ou une preuve de ce fait.

D'accord, cela devrait être suffisant explication de fond. Retour à votre (modifié) après. Vous dites « Le erf (z) scipy.special me demander de définir exactement ce que t est initialement ». Je ne sais pas ce que vous entendez par là; où est-<=> (temps?) entrer dans ce tout? Espérons que l'explication a démystifié au-dessus de la fonction d'erreur un peu et il est maintenant plus clair que la raison pour laquelle la fonction d'erreur est la fonction pour le travail.

Votre code Python est OK, mais je préfère une fermeture sur une 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

L'utilisation d'une fermeture permet précalcul des constantes et <=> <=>, de sorte que la fonction de retour devra faire moins de travail à chaque fois qu'il est appelé (qui peut être important si vous intégrez, ce qui signifie qu'il sera appelé à plusieurs reprises). J'ai également évité toute utilisation de l'opérateur d'exponentiation <=>, qui est plus lente que l'écriture seulement la mise au carré en panne et hissé la fracture de la boucle interne et l'a remplacée par une multiplication. Je ne l'ai pas regardé tout à leur mise en œuvre en Python, mais de mon dernier réglage d'une boucle interne pour la vitesse pure en utilisant l'assemblage x87 cru, il me semblede se rappeler que additions, soustractions, ou multiplications prendre environ 4 cycles CPU chacun, divise environ 36 et exponentiation environ 200. C'était il y a quelques années, alors prenez ces chiffres avec un grain de sel; encore, il illustre leur complexité relative. De plus, le calcul de la manière-<=> force brute est une très mauvaise idée; il y a des trucs que vous pouvez prendre lors de l'écriture d'une bonne mise en œuvre de qui le rendent <=> nettement plus rapide et plus précis qu'un exponentiation style général <=>.

Je ne l'ai jamais utilisé la version numpy des constantes pi et e; Je me suis toujours coincé avec les versions du module mathématique ancienne plaine. Je ne sais pas pourquoi vous préférez peut-être un ou l'autre.

Je ne sais pas ce que vous allez pour l'appel <=>. Devrait intégrer <=> une gaussienne renormalisée de moins l'infini à plus l'infini, et doit toujours cracher 10 (votre facteur de normalisation), puisque la gaussienne intègre à 1 sur la ligne réelle. Toute réponse loin de 10 (je n'attendre exactement 10 puisque seulement une <=> approximation, après tout) signifie quelque chose est foiré quelque part ... difficile de dire ce qui se visse sans connaître la la valeur de retour réelle et, éventuellement, le fonctionnement interne de <=>.

Il faut espérer que a démystifié la confusion, et expliqué pourquoi la fonction d'erreur est la bonne réponse à votre problème, ainsi que la façon de le faire vous-même si vous êtes curieux. Si l'une de mon explication n'a pas été claire, je vous suggère de jeter un oeil rapide sur Wikipédia premier; si vous avez encore des questions, ne hésitez pas à demander.

Answer 2

navires SciPy avec la fonction "d'erreur", intégrale gaussienne aka:

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

Answer 3

Je suppose que vous manipulez gaussiennes à plusieurs variables; le cas échéant, SciPy a déjà la fonction que vous cherchez: il est appelé MVNDIST ( "multivariée de distribution normale) La documentation SciPy est, comme toujours, terrible, donc je ne peux pas trouver, même si la fonction est enterré, mais il est là quelque part . la documentation est facilement le pire des SciPy, et m'a frustré sans fin dans le passé.

simple variable gaussiennes il suffit d'utiliser la bonne vieille fonction d'erreur, dont de nombreuses implémentations sont disponibles.

Quant à attaquer le problème en général, oui, comme James Thompson mentionne, vous voulez juste écrire votre propre fonction de distribution gaussienne et le nourrir à quatre (). Si vous pouvez éviter l'intégration généralisée, cependant, il est une bonne idée de le faire - des techniques d'intégration spécialisés pour une fonction particulière (comme MVNDIST utilise) vont être beaucoup plus rapide que l'intégration multidimensionnelle Monte Carlo standard, qui peut être extrêmement lent pour une grande précision.

Answer 4

La distribution gaussienne est aussi appelée une distribution normale. La fonction cdf dans le module de norme scipy fait ce que vous voulez.

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

Pourquoi savez pas seulement toujours faire votre intégration de -infinity à + l'infini, de sorte que vous avez toujours la réponse? (Blague!)

Je pense que la seule raison qu'il n'y a pas déjà une gaussienne en conserve dans SciPy est qu'il est une fonction triviale à écrire. Votre suggestion d'écrire votre propre fonction et le transmettre à quad d'intégrer des sons excellents. Il utilise l'outil SciPy accepté pour ce faire, il est l'effort de code minimal pour vous, et il est très facile à lire pour les autres, même si elles ont jamais vu SciPy.

Que voulez-vous dire par un intégrateur de largeur fixe? Voulez-vous dire en utilisant un algorithme différent de ce que QUADPACK utilise?

Edit: Pour être complet, voici quelque chose comme ce que je vais essayer d'une gaussienne avec la moyenne de 0 et de 1 écart-type de 0 à + infini:

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 )

C'est un peu laid parce que la fonction de Gauss est un peu long, mais encore assez trivial à écrire.