Verhindern von Unterschreitungen als Logarithmus der Wahrscheinlichkeitsberechnungs, die einer normalen Probe in Python in einem bestimmten Intervall fällt

Question

Gegeben seien zwei reellen Zahlen x und y, möchte ich die folgende Funktion in Python berechnen:

log Pr [ x <= t <= y ], 

wo t abgetastet wird aus einer normalen Verteilung.

Eine naive Implementierung ist scipy.stats.norm zu verwenden.

np.log(scipy.stats.norm.cdf(y) - scipy.stats.norm.cdf(x)) 

Leider ist dies führt zu einer Unterschreitung, wenn x und y weit von 0 sind. Wie kann man einen solchen numerischen Fehler verhindern?

Answer 1

Dieses Problem ist viel stabiler, wenn im Logspace erfolgt.

Der Trick besteht darin, scipy.stats.norm.logcdf für Werte kleiner als Null und scipy.stats.norm.logsf für Werte größer als Null zu verwenden.

Diese mit einem stabilen Algorithmus kombiniert log(exp(y) - exp(x)) vernünftig gibt Ergebnisse für die Berechnung

import numpy as np 
from scipy.stats import norm 

def log_subtract(x, y): 
    return x + np.log1p(-np.exp(y-x)) 

def lnprob(x, y): 
    if x < 0: 
     return log_subtract(norm.logcdf(y), norm.logcdf(x)) 
    else: 
     return log_subtract(norm.logsf(x), norm.logsf(y))