欺骗 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) uses a Taylor series expansion of the error function to iteratively compute the log CDF directly，这比简单取log(ndtr(x)).

数值稳定得多

更新

正如您在评论中提到的，如果 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

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

您需要计算, which you correctly identified would be 的积分。打开括号您可以将求和的两个部分合并。

Scipy 有这个 imaginary error function implemented

第二部分更难：

这是 generalized hypergeometric function. Sadly it looks like scipy does not have an implementation of it, but this package 声称的。

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

Answer 3

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

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

Tricking numpy/python into representing very large and very small numbers

python

floating-point

numbers

numpy

scipy

更新