欺骗numpy/python代表非常大和非常小的数字

Question

我需要在从低到低的范围内计算以下函数的积分-150:

import numpy as np
from scipy.special import ndtr

def my_func(x):
    return np.exp(x ** 2) * 2 * ndtr(x * np.sqrt(2))

问题是这部分功能

np.exp(x ** 2)

倾向于无穷大 - 我得到inf的值x小于约-26.

而这部分功能

2 * ndtr(x * np.sqrt(2))

这相当于

from scipy.special import erf

1 + erf(x)

倾向于0.

因此,一个非常非常小的数字非常非常大,应该给我一个合理大小的数字 - 但是,而不是那个,python给了我nan.

我能做些什么来规避这个问题？

Answer 1

我认为@ askewchan的解决方案的组合scipy.special.log_ndtr将会做到这一点:

from scipy.special import log_ndtr

_log2 = np.log(2)
_sqrt2 = np.sqrt(2)

def my_func(x):
    return np.exp(x ** 2) * 2 * ndtr(x * np.sqrt(2))

def my_func2(x):
    return np.exp(x * x + _log2 + log_ndtr(x * _sqrt2))

print(my_func(-150))
# nan

print(my_func2(-150)
# 0.0037611803122451198

因为x <= -20,log_ndtr(x) 使用误差函数的泰勒级数展开来直接迭代地计算log CDF,这在数值上比简单地更稳定log(ndtr(x)).

更新

正如您在评论中提到的,exp如果x足够大,也会溢出.虽然你可以解决这个问题mpmath.exp,但是一个更简单,更快速的方法就是np.longdouble在我的机器上投射到最多1.189731495357231765e + 4932的值:

import mpmath

def my_func3(x):
    return mpmath.exp(x * x + _log2 + log_ndtr(x * _sqrt2))

def my_func4(x):
    return np.exp(np.float128(x * x + _log2 + log_ndtr(x * _sqrt2)))

print(my_func2(50))
# inf

print(my_func3(50))
# mpf('1.0895188633566085e+1086')

print(my_func4(50))
# 1.0895188633566084842e+1086

%timeit my_func3(50)
# The slowest run took 8.01 times longer than the fastest. This could mean that
# an intermediate result is being cached  100000 loops, best of 3: 15.5 µs per
# loop

%timeit my_func4(50)
# The slowest run took 11.11 times longer than the fastest. This could mean
# that an intermediate result is being cached  100000 loops, best of 3: 2.9 µs
# per loop

Answer 2

已经有这样一个函数：erfcx。我认为erfcx(-x)应该给你你想要的被积函数（注意1+erf(x)=erfc(-x)）。