Der beste Weg, eine Python-Funktion zu schreiben, die eine Gaußsche integriert?

StackOverflow https://stackoverflow.com/questions/509994

Frage

Bei dem Versuch, scipy Quad-Methode eine Gaußsche zu integrieren (können sagen, es gibt eine Gaußsche Methode namens Gauss), ich hatte Probleme benötigten Parameter übergeben und Quad die Integration über die richtige Variable zu tun verlassen gauss. Hat jemand ein gutes Beispiel davon, wie Quad verwenden w / eine mehrdimensionale Funktion?

Aber das führte mich zu einer großen Frage nach dem besten Weg, eine Gaußsche im Allgemeinen zu integrieren. Ich habe keine Gaußsche integrieren in scipy (zu meiner Überraschung) finden. Mein Plan war, eine einfache Gaußsche Funktion zu schreiben, und übergeben Sie es (oder vielleicht jetzt eine feste Breite Integrator) auf Quad. Was würden Sie tun?

Edit:. Mit fester Breite etwas wie trapz Bedeutung, die eine feste dx verwendet Bereiche unter einer Kurve zu berechnen

Was ich bin gekommen, um so weit ist ein Verfahren 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 machen und Varianz Ich muss vor der Integration. 

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 genannt werden muss) und 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 nicht notwendigerweise übereinstimmen N = 10, sigma = 2, mu = 0, die den mehr erweiterte Definition aufgefordert werden.

Die erf (z) in scipy.special würde mich erfordern, genau zu definieren, was t ist zunächst, aber es ist schön zu wissen, dass es da ist.

War es hilfreich?

Lösung

Okay, scheinen Sie ziemlich verwirrt über einige Dinge zu sein. Lassen Sie uns am Anfang beginnen: Sie haben eine „mehrdimensionale Funktion“ erwähnt, aber dann gehen Sie auf die übliche variable Gauß-Kurve zu diskutieren. Dies ist nicht eine mehrdimensionale Funktion: Wenn Sie es integrieren, können Sie nur eine Variable integrieren (x). Diese Unterscheidung ist wichtig zu machen, weil es ist ein Monster genannt „multivariate Gauß-Verteilung“, das eine wahre multidimensionale Funktion ist und, falls integriert, erfordert mehr als zwei oder mehr Variablen zu integrieren (die die teueren Monte verwendet Carlo Technik, die ich bereits erwähnt). Aber Sie scheinen nur über die regulären variablen Gaußschen sprechen, die viel einfacher zu handhaben ist, zu integrieren, und so weiter.

Die eine variable Gauß-Verteilung hat zwei Parameter, sigma und mu und ist eine Funktion einer einzigen Variable, die wir x bezeichnen werde. Sie erscheinen auch um einen Normalisierungsparameter n tragen werden (die in einigen Anwendungen nützlich ist). Normalisierungsparameter sind in der Regel nicht , die in Berechnungen, da man sie nur am Ende wieder auf heften kann (denken Sie daran, die Integration ist ein linearer Operator: int(n*f(x), x) = n*int(f(x), x)). Aber wir können es sich herumtragen, wenn Sie mögen; 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))

(lesen, dass "die Normalverteilung von x gegeben sigma, mu und n ist gegeben durch ...") So weit, so gut; dies entspricht der Funktion, die Sie haben. Beachten Sie, dass der einzige true Variable hier ist x. Die anderen drei Parameter sind fixed für eine bestimmte Gaussian

Jetzt für eine mathematische Tatsache: Es ist beweisbar wahr, dass alle Gauß-Kurven, die die gleiche Form haben, sind sie um ein wenig nur verschoben. So können wir mit N(x|0,1,1) arbeiten, die „Standardnormalverteilung“ genannt, und nur unsere Ergebnisse zurück an die allgemeinen Gauß-Kurve übersetzen. Wenn Sie also das Integral von N(x|0,1,1) haben, können Sie trivialer das Integral jeden Gaußschen berechnen. Dieses Integral erscheint so häufig, dass es einen besonderen Namen hat: die Fehlerfunktion erf. Aufgrund einiger alten Konventionen, es ist nicht genau erf; es gibt ein paar additiven und multiplikativen Faktoren auch um getragen zu werden.

Wenn Phi(z) = integral(N(x|0,1,1), -inf, z); das heißt, Phi(z) das Integral der Standardnormalverteilung von minus unendlich bis z, dann ist es durch die Definition der Fehlerfunktion, dass

true

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

Auch 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 mu ist, sigma und n von minus unendlich bis z, dann ist es wahr, durch die Definition der Fehlerfunktion, die

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

Hier finden Sie aktuelle , wenn Sie mehr Detail wollen oder ein Beweis für diese Tatsache.

Okay, das sollte genug Hintergrund Erklärung. Zurück zu Ihrem (bearbeitet) Post. Sie sagen: „Die erf (z) in scipy.special würde mich erfordern, genau zu definieren, was t zunächst ist“. Ich habe keine Ahnung, was Sie damit meinen; wo kommt t (Zeit) überhaupt in diese eingeben? Hoffentlich wird die Erklärung über entmystifiziert hat die Fehlerfunktion ein wenig und es ist jetzt klarer, warum die Fehlerfunktion ist die richtige Funktion für den Job.

Ihr Python-Code ist in Ordnung, aber ich würde einen Verschluss über eine 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 der zurück Funktion wilIch brauche weniger Arbeit jedes Mal, es zu tun haben genannt (was wichtig sein kann, wenn Sie es sind zu integrieren, was bedeutet es oft genannt werde). Außerdem habe ich jede Verwendung des Exponential-Operator ** vermieden, was als langsamer ist nur das Schreiben der Quadrierung heraus und gehißt divide die aus der inneren Schleife und ersetzen sie durch eine mehrfach. Ich habe überhaupt nicht auf ihrer Umsetzung in Python gesucht, aber von meiner letzten Zeit eine innere Schleife für reine Geschwindigkeit unter Verwendung von roher x87 Montage Tuning, ich mich zu erinnern, dass addiert, subtrahiert, oder vervielfachen nimmt etwa 4 CPU-Zyklen jeweils teilt über 36 und Potenzierung etwa 200. Das vor ein paar Jahren war, so die Zahlen mit einem Körnchen Salz; Dennoch zeigt sie ihre relative Komplexität. Wie gut, ist die Brute-Force-Weg exp(x) Berechnung eine sehr schlechte Idee; es gibt Tricks, die Sie nehmen können, wenn eine gute Umsetzung von exp(x) zu schreiben, die es als eine allgemeine a**b Stil Potenzierung wesentlich schneller und genauer zu machen.

Ich habe noch nie verwendet die numpy Version der Konstanten pi und e; Ich habe immer mit der einfachen alten Mathe-Moduls Versionen stecken. Ich weiß nicht, warum Sie entweder ein vorziehen könnten.

Ich bin nicht sicher, was Sie mit dem quad() Anruf gehen. quad(gen_gauss, -inf, inf, (10,2,0)) sollte eine renormierte Gaussian von minus unendlich bis plus unendlich integrieren, und soll immer ausspucken 10 (Ihren Normierungsfaktor), da die Gauß-1 über die durchgezogene Linie integriert. Jede Antwort weit von 10 (Ich würde erwarten, nicht genau 10 seit quad() ist nur eine Annäherung, nachdem alle) bedeutet, ist etwas irgendwo vermasselt ... schwer zu sagen, was nach oben geschraubt ist, ohne die tatsächliche Rückkehr zu wissen Wert und möglicherweise das Innenleben von quad().

Hoffentlich hat einige der Verwirrung entmystifiziert, und erklärt, warum die Fehlerfunktion die richtige Antwort auf Ihr Problem ist, und wie es alles selbst tun, wenn Sie neugierig sind. Wenn einer meiner Erklärung nicht klar war, schlage ich zuerst einen kurzen Blick auf Wikipedia nehmen; Wenn Sie noch Fragen haben, zögern Sie nicht zu fragen.

Andere Tipps

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

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

Ich nehme an, Sie multivariate Gaussians sind Handhabung; wenn ja, hat SciPy bereits die Funktion Sie suchen: es heißt MVNDIST ( "Multivariate Normalverteilung) Die SciPy Dokumentation ist, wie immer, schrecklich, so kann ich nicht einmal, wo die Funktion begraben, aber es ist dort irgendwo . die Dokumentation ist einfach der schlimmste Teil SciPy und hat mich in der Vergangenheit zu keinem Ende frustriert.

Single-Variable Gaussians benutzen Sie einfach die gute alte Fehlerfunktion, von denen verfügbar viele Implementierungen sind.

Wie für das Problem im allgemeinen Angriff, ja, wie James Thompson erwähnt, möchten Sie gerade Ihre eigene Gaußschen Verteilungsfunktion schreiben und ihn an quad (). Wenn Sie die verallgemeinerte Integration vermeiden können, aber es ist eine gute Idee, dies zu tun - spezialisierte Integrationstechniken für eine bestimmte Funktion (wie MVNDIST verwendet) werden viel schneller als eine Standard-Monte Carlo mehrdimensionale Integration, die extrem langsam sein kann für hohe Genauigkeit.

Die Gaußsche Verteilung ist auch eine Normalverteilung genannt. Die CDF-Funktion in der scipy Norm Modul tut, 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

Warum nicht einfach immer tun, um Ihre Integration von -unendlich bis + unendlich, so dass Sie immer wissen, die Antwort? (Scherz!)

Meine Vermutung ist, dass der einzige Grund, das ist es nicht bereits eine Dose Gauß-Funktion in SciPy ist, dass es eine triviale Funktion zu schreiben. Ihr Vorschlag über eine eigene Funktion zu schreiben und weitergeben klingt hervorragend Quad zu integrieren. Es nutzt das akzeptierte SciPy Werkzeug, dies zu tun, ist es minimal Code Aufwand für Sie, und es ist sehr gut lesbar für andere Menschen, auch wenn sie noch nie SciPy gesehen.

Was genau meinst du mit einem Fest Breiten Integrator? Meinen Sie einen anderen Algorithmus als was auch immer QUADPACK verwendet?

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 )

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

Lizenziert unter: CC-BY-SA mit Zuschreibung
Nicht verbunden mit StackOverflow
scroll top