Оценить журнал (1 - normal_cdf (x)) численно стабильным способом

Question

Как оценить log(1 - normal_cdf(x)) численно устойчивым способом? Здесь normal_cdf - кумулятивная функция распределения стандартного нормального распределения.

Например, в Python:

import scipy 
from scipy.stats import norm

np.log(1 - norm.cdf(10))

дает -inf с RuntimeWarning: divide by zero encountered in log, поскольку norm.cdf(10) почти равно 1. Есть ли такая функция, как logsumexp, которая может избежать числового занижения?

Answer 1

@ Предложение HongOoi воспользоваться симметрией велико. Но для произвольного распределения в scipy.stats (включая norm) вы можете использовать метод logsf именно для этого вычисления. sf обозначает функция выживания , которая является названием функции 1 - cdf(x).

Например,

In [25]: import numpy as np

In [26]: from scipy.stats import norm, gamma

Вот пример norm.logsf:

In [27]: norm.logsf(3, loc=1, scale=1.5)
Out[27]: -2.3945773661586434

In [28]: np.log(1 - norm.cdf(3, loc=1, scale=1.5))
Out[28]: -2.3945773661586434

А вот пример gamma.logsf:

In [29]: gamma.logsf(1.2345, a=2, scale=1.8)
Out[29]: -0.16357333194167956

In [30]: np.log(1 - gamma.cdf(1.2345, a=2, scale=1.8))
Out[30]: -0.16357333194167956

Это показывает, почему нужно использовать logsf(x) вместо log(1 - cdf(x)):

In [35]: norm.logsf(50, loc=1, scale=1.5)
Out[35]: -537.96178420294677

In [36]: np.log(1 - norm.cdf(50, loc=1, scale=1.5))
/Users/warren/miniconda3scipy/bin/ipython:1: RuntimeWarning: divide by zero encountered in log
  #!/Users/warren/miniconda3scipy/bin/python
Out[36]: -inf

Answer 2

Поскольку нормальное распределение симметрично относительно 0, мы имеем

1 - F(x) = P(X > x)
         = P(X < -x)
         = F(-x)

Следовательно

np.log(1 - norm.cdf(10)) = np.log(norm.cdf(-10))
                         = norm.logcdf(-10)