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

趋向于无穷大 -- 对于 x 的值小于大约 -26，我得到 inf。

还有这部分功能

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)uses a Taylor series expansion of the error function to iteratively compute the log CDF directly，这比简单地采用log(ndtr(x))在数值上稳定得多。

---

更新

正如您在 cmets 中提到的，如果 x 足够大，exp 也会溢出。虽然您可以使用 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)）。

Answer 3

不确定这会有多大帮助，但这里有一些想法太长了，无法发表评论。

您需要计算 的积分，您的correctly identified 将是。打开括号，您可以整合总和的两个部分。

Scipy 有这个imaginary error function implemented

第二部分更难：

这是generalized hypergeometric function。可悲的是，它看起来像 scipy does not have an implementation of it，但 this package 声称确实如此。

这里我使用了不带常数的不定积分，知道了from to 的值，就很清楚如何使用定积分了。