Wenn ich eine Zufallszahl z habe, die als Summe von zwei anderen Zufallsnummern definiert ist, x und y , dann ist die Wahrscheinlichkeitsverteilung von z die Faltung der Wahrscheinlichkeitsverteilungen für x und y . Die Faltung ist im Grunde auf ein Integral des Produkts der Vertriebsfunktionen. Oft gibt es keine analytische Lösung für das Integral in der Faltung, daher muss sie mit einem grundlegenden Quadraturalgorithmus berechnet werden. In Pseudocode: < /p>
Code: Select all
prob_z(z) = integrate(lambda t: prob_x(t) * prob_y(z-t), -inf, inf)
Für ein konkretes Beispiel kann die Summe z einer normal verteilten Variablen x und ein normalerweise verteiltes Protokollvariable y mit dem folgenden Python/Scipy -Code berechnet werden:
Code: Select all
from scipy.integrate import quad
from scipy.stats import norm, lognorm
from scipy import log
prob_x = lambda x: norm.pdf(x, 0, 1) # N(mu=0, sigma=1)
prob_y = lambda y: lognorm.pdf(y, 0.1, scale=10) # LogN(mu=log(10), sigma=0.1)
def prob_z(z):
return quad(lambda t: prob_x(t)*prob_y(z-t), -inf, inf)
< /code>
Jetzt möchte ich die Protokollwahrscheinlichkeit berechnen. Die naive Lösung ist einfach zu tun: < /p>
def log_prob_z(z):
return log(prob_z(z))
< /code>
Dies ist jedoch numerisch instabil. Nach etwa 39 Standardabweichungen betragen die Wahrscheinlichkeitsverteilungen numerisch 0,0. Selbst wenn die Protokollwahrscheinlichkeit einen endlichen Wert aufweist, kann sie nicht außerhalb dessen berechnet werden, indem einfach das Protokoll der Wahrscheinlichkeit eingenommen wird. Vergleiche norm.pdf (39, 1, 0)
, die 0,0 bis norm.Logpdf (39, 1, 0) ist, was ungefähr -761 ist. Offensichtlich berechnet Scipy logpdf als log (pdf) -es findet einen anderen Weg -da es sonst -inf eine minderwertige Antwort zurückgibt. Ebenso möchte ich einen anderen Weg für mein
Problem finden. /> Die Frage lautet: Weiß jemand, wie ich Protokoll (Quad (...)) neu ordnen kann.>
Mobile version