Предотвращения недорасхода при вычислении логарифма вероятности того, что нормальный образец попадет в определенном интервале в питоне

Question

Учитывая два реала x и y, я хочу, чтобы вычислить следующую функцию в Python:

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

где t дискретизируются из нормального распределения.

Одна наивная реализация заключается в использовании scipy.stats.norm.

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

К сожалению, это вызывает опустошение, когда x и y далеки от 0. Как предотвратить такую ​​численную ошибку?

Answer 1

Эта проблема намного более стабильна, если она выполнена в logspace.

Хитрость заключается в использовании scipy.stats.norm.logcdf для значений, меньших нуля, и scipy.stats.norm.logsf для значений, больших нуля.

Это в сочетании со стабильным алгоритмом расчета log(exp(y) - exp(x)) дает приемлемые результаты

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))