Der beste Weg, um eine Python-Funktion zu schreiben, die einen Gaussian integriert?

Question

Bei dem Versuch, die quad-Methode von scipy zu verwenden, um einen Gaussian zu integrieren (ich nehme an, es gibt eine Gauss-Methode namens Gauss), hatte ich Probleme, benötigte Parameter an Gauss zu übergeben und Quad zu verlassen, um die Integration über die korrekte Variable durchzuführen. Hat jemand ein gutes Beispiel dafür, wie man Quad mit einer multidimensionalen Funktion benutzt?

Aber das führte mich zu einer größeren Frage über den besten Weg, eine Gaussian im Allgemeinen zu integrieren. Ich habe kein Gaussian in scipy integriert (zu meiner Überraschung). Mein Plan war, eine einfache Gauss-Funktion zu schreiben und sie an Quad (oder vielleicht einen Integrator mit fester Breite) weiterzuleiten. Was würden Sie tun?

Bearbeiten: Feste Breite bedeutet etwas wie Trapz, die eine feste DX verwendet, um Bereiche unter einer Kurve zu berechnen.

Was ich bis jetzt gekommen bin, ist eine Methode make___gauss, die eine Lambda-Funktion zurückgibt, die dann in Quad gehen kann. Auf diese Weise kann ich eine normale Funktion mit dem Durchschnitt und der Varianz machen, die ich vor der Integration benötige.

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) 

Als ich versuchte, eine allgemeine Gaußsche Funktion übergeben (die mit x, N, mu, und Sigma aufgerufen werden muss), und das Füllen in einige der Werte unter Verwendung von Quad wie

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

die Parameter 10, 2 und 0 stimmen NICHT notwendigerweise mit N = 10, sigma = 2, mu = 0 überein, was die erweiterte Definition anregte.

Das erf (z) in scipy.special würde erfordern, dass ich genau definiere, was t ist, aber es ist schön zu wissen, dass es da ist.

Answer 1

Okay, scheinen Sie ziemlich verwirrt über einige Dinge zu sein. Fangen wir gleich am Anfang an: Sie haben eine "multidimensionale Funktion" erwähnt, diskutieren dann aber die übliche Gauss-Kurve mit einer Variablen. Dies ist nicht eine multidimensionale Funktion: Wenn Sie es integrieren, integrieren Sie nur eine Variable (x). Die Unterscheidung ist wichtig zu machen, weil ein Monster ist, das eine "multivariate Gaußsche Verteilung" genannt wird, die eine echte multidimensionale Funktion ist und, wenn integriert, Integration über zwei oder mehr Variablen erfordert (die die teure Monte-Carlo-Technik verwendet, die ich zuvor erwähnt habe)). Aber Sie scheinen nur über die reguläre Gaußsche Variable mit einer Variablen zu sprechen, die viel einfacher zu bearbeiten und zu integrieren ist.

Die Gauss-Verteilung mit einer Variablen hat zwei Parameter, sigma und mu, und ist eine Funktion einer einzelnen Variablen, die wir x bezeichnen. Sie scheinen auch einen Normalisierungsparameter n herumzutragen (was in einigen Anwendungen nützlich ist). Normalisierungsparameter sind normalerweise nicht in den Berechnungen enthalten, da Sie sie am Ende nur wieder anheften können (denken Sie daran, Integration ist ein linearer Operator: int(n*f(x), x) = n*int(f(x), x)). Aber wir können es herumtragen, wenn Sie möchten; die Notation ich für eine Normalverteilung wie ist dann

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

(das als gelesen „die Normalverteilung von xsigma, mu gegeben, und n ist gegeben durch ...“) So weit, so gut; Dies entspricht der Funktion, die Sie haben. Beachten Sie, dass die einzige wahre Variable ist hier x: die anderen drei Parameter sind festen für eine bestimmte Gaußsche.

Jetzt für eine mathematische Tatsache: es ist nachweislich wahr, dass alle Gaußkurven die gleiche Form haben, sie sind nur ein wenig verschoben. So können wir mit N(x|0,1,1) arbeiten, die "Standardnormalverteilung" genannt wird, und unsere Ergebnisse einfach in die allgemeine Gaußkurve zurückübersetzen. Wenn Sie also das Integral N(x|0,1,1) haben, können Sie das Integral beliebiger Gaussiane trivial berechnen. Dieses Integral erscheint so häufig, dass es einen speziellen Namen hat: die Fehlerfunktionerf. Wegen einiger alter Konventionen ist es nicht genauerf; es gibt ein paar additive und multiplikative Faktoren, die auch herumgetragen werden.

Wenn ; das heißt, Phi(z) das Integral der Normalverteilung von minus unendlich bis zu z, dann ist es durch die Definition der Fehlerfunktion wahr, dass

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

Ebenso, wenn Phi(z | mu, sigma, n) = integral(N(x|sigma, mu, n), -inf, z); das heißt, Phi(z | mu, sigma, n) das Integral der Normalverteilung gegebenen Parameter ist mu, sigma und n von minus unendlich bis z, dann ist es durch die Definition der Fehlerfunktion, dass

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

Werfen Sie einen Blick auf the Wikipedia article on the normal CDF, wenn Sie mehr Details oder einen Beweis für diese Tatsache wollen.

Okay, das sollte genug Hintergrunderklärung sein. Zurück zu Ihrem (redigierten) Beitrag. Sie sagen "Das erf (z) in scipy.special würde erfordern, dass ich genau definiere, was t anfänglich ist". Ich habe keine Ahnung, was du damit meinst; Wohin geht t (Zeit?) überhaupt ein? Hoffentlich hat die obige Erklärung die Fehlerfunktion etwas entmystifiziert und es ist jetzt klarer, warum die Fehlerfunktion die richtige Funktion für den Job ist.

Ihr Python-Code ist in Ordnung, aber ich würde einen Verschluss über einen Lambda bevorzugen:

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 

einen Verschluss Mit ermöglicht precomputation von Konstanten k und s, so dass die zurück Funktion benötigt weniger tun jedes Mal arbeiten es heißt (was kann wichtig sein, wenn Sie es integrieren, was bedeutet, dass es viele Male aufgerufen wird). Außerdem habe ich jede Verwendung des Potenzierungsoperators ** vermieden, der langsamer ist als nur das Quadrieren zu schreiben, und die Teilung aus der inneren Schleife herausgehoben und durch eine Multiplikation ersetzt. Ich habe noch nicht ihre Implementierung in Python betrachtet, aber seit ich das letzte Mal eine innere Schleife für reine Geschwindigkeit unter Verwendung roher x87-Assembly eingestellt habe, denke ich, dass Addieren, Subtrahieren oder Multiplizieren jeweils etwa 4 CPU-Zyklen benötigen 36, und Potenzierung etwa 200. Das war vor ein paar Jahren, also nimm diese Zahlen mit einem Körnchen Salz; Dennoch zeigt es ihre relative Komplexität. Auch die Berechnung der exp(x) Brute-Force-Methode ist eine sehr schlechte Idee; Es gibt Tricks, die Sie beim Schreiben einer guten Implementierung von exp(x), die es deutlich schneller und genauer als eine allgemeine Exponentiation a**b machen können.

Ich habe nie die numpy Version der Konstanten pi und e verwendet; Ich habe mich immer an die einfachen Versionen des alten Mathematikmoduls gehalten. Ich weiß nicht, warum du eines bevorzugen solltest.

Ich bin mir nicht sicher, was Sie mit dem quad() Aufruf gehen. quad(gen_gauss, -inf, inf, (10,2,0)) sollte ein renormiertes Gauß-Signal von minus unendlich bis plus unendlich integrieren und sollte immer 10 (Ihren Normalisierungsfaktor) ausspucken, da sich die Gaußsche Kurve über die reale Linie in 1 integriert. Jede Antwort weit weg von 10 (ich würde nicht erwarten genau 10 seit quad() ist nur eine Annäherung, schließlich) bedeutet etwas ist irgendwo vermasselt ... schwer zu sagen, was vermasselt ist, ohne den tatsächlichen Rückgabewert und möglicherweise die innere wissen Arbeiten von quad().

Hoffentlich hat das einige der Verwirrung entmystifiziert und erklärt, warum die Fehlerfunktion die richtige Antwort auf Ihr Problem ist, und wie Sie alles selbst machen können, wenn Sie neugierig sind. Wenn eine meiner Erklärungen nicht klar war, schlage ich vor, zuerst einen kurzen Blick auf Wikipedia zu werfen; Wenn Sie noch Fragen haben, zögern Sie nicht zu fragen.

Answer 2

Warum sollten Sie nicht immer Ihre Integration von -unendlich bis + unendlich machen, damit Sie immer die Antwort wissen? (Scherz!)

Meine Vermutung ist, dass der einzige Grund, dass es nicht bereits eine vordefinierte Gaußsche Funktion in SciPy gibt, dass es eine triviale Funktion ist, zu schreiben. Ihr Vorschlag über das Schreiben Ihrer eigenen Funktion und das Übergeben an Quad, um Sounds zu integrieren, ist ausgezeichnet. Es verwendet das akzeptierte SciPy-Tool, um dies zu tun, es ist ein minimaler Code-Aufwand für Sie, und es ist sehr lesbar für andere Menschen, auch wenn sie SciPy noch nie gesehen haben.

Was genau meinen Sie mit einem Integrator mit fester Breite? Meinst du einen anderen Algorithmus als den, den QUADPACK benutzt?

Edit: Für Vollständigkeit, hier ist etwas, was ich für eine Gaußsche mit dem Mittelwert von 0 und eine Standardabweichung von 1 von 0 bis + unendlich versuchen würde:

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) 

, dass, weil der Gauß ein wenig hässlich Funktion ist ein wenig lang, aber immer noch ziemlich trivial zu schreiben.

Answer 3

Ich gehe davon aus, dass Sie mit multivariaten Gaussians umgehen; Wenn ja, SciPy hat bereits die Funktion, nach der Sie suchen: Sie heißt MVNDIST ("MultiVariate Normal Distribution"). Die SciPy-Dokumentation ist wie immer schrecklich, sodass ich nicht einmal herausfinden kann, wo die Funktion verborgen ist, sondern it's in there somewhere. die Dokumentation ist einfach der schlimmste Teil SciPy und hat mich in der Vergangenheit zu keinem Ende frustriert.

Einzel-variablen Gaussians benutzen Sie einfach die gute alte Fehlerfunktion, von denen viele Implementierungen zur Verfügung stehen.

Was das Problem allgemein anzugreifen, ja, wie James Thompson erwähnt, wollen Sie einfach Ihre eigene Gauss'sche Verteilungsfunktion schreiben und sie an quad() übergeben.Wenn Sie jedoch die verallgemeinerte Integration vermeiden können, ist dies eine gute Idee - spezialisierte Integrationstechniken für eine bestimmte Funktion (wie MVNDIST verwendet) werden viel schneller sein als eine standardmäßige Monte-Carlo-multidimensionale Integration, die extrem langsam sein kann für hohe Genauigkeit.

Answer 4

scipy Schiffe mit der "Fehlerfunktion", auch bekannt als Gaußsche Integral:

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

Answer 5

Die Gaußverteilung wird auch als Normalverteilung bezeichnet. Die cdf-Funktion im scipy-Norm-Modul erfüllt das, was Sie wollen.

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