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

Question

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

9

我需要计算以下函数的积分，范围从-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。

我应该怎么做来避免这个问题？

- abcd

@ReblochonMasque 不是，你知道我可以在哪里找到吗？我肯定没有数学技能自己解决它。 - abcd

可能是mathstackexchange - 或者wolframalpha - 或者sympy。 - Reblochon Masque

1

这样的东西有帮助吗？np.exp(x**2 + np.log(2) + np.log(ndtr(x*np.sqrt(2)))) - askewchan

1

这是否会实质性地影响积分的值？或者，您可以通过分析计算出log(ndtr(x))是什么，并像 exp 项一样进行拆分...您明白了。 - askewchan

1

你可以使用scipy.special.log_ndtr来消除@askewchan解决方案中的log调用。 - ali_m

显示剩余6条评论

3个回答

4

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

- Simon Byrne

erf(x)是一种奇函数：erf(-x) = -erf(x)。因此，erfc(-x) = 1-erf(-x) = 1 + erf(x)（第一个等式是erfc的定义，第二个等式使用了奇对称性）。 - Warren Weckesser

请注意，通过翻转参数的符号，您可以将quad(my_func, -14, -4)替换为quad(erfcx, 4, 14)。 - Warren Weckesser

@WarrenWeckesser 谢谢。最后一个问题：erfcx 定义为 exp(x ** 2) * erfc(x)。如果我执行 erfcx(-x)，那么我不是得到了 exp(-x ** 2) * erfc(-x)，而我想要的是 exp(x ** 2) * erfc(-x) 吗？或者，等等，我想它会给我 exp((-x) ** 2) * erf((-x))，这正是我想要的。 - abcd

2

“不确定这对你有多大帮助，但这里有几个想法，它们太长了无法在评论中说出。

您需要计算 $2 \cdot e^{x^2} \cdot f(\sqrt{2}x)$ 的积分，您已经正确地确认它将会是 $e^{x^2}*(1 + erf(x))$ 。展开括号后，您可以积分求解两部分的和。”

Scipy已经实现了这个虚数误差函数。

第二部分比较困难：

这是一个广义超几何函数。不幸的是，看起来scipy 没有实现它，但这个软件包声称有。

在这里，我使用了没有常数的不定积分，知道了from to的值，就可以清楚地如何使用定积分。

- Salvador Dali

似乎困扰我这种策略的是 scipy.special.erfi(50) 的计算结果为 inf。 - abcd

数学问题：如果我的积分是“半确定”的，即我有一个上限，但我的下限是-inf，我该如何处理？我不确定如何在-inf处计算您问题中的函数。 - abcd

@dbliss scipy.special.erfi(50) 的结果为无穷大。如果b=超级大数值，那么很可能b + a的结果也是超级大的数，其中a>0（因为当x>1时，第二个积分也为正）。不太确定您想如何计算它。 - Salvador Dali

(1) 是的，但是 mpmath.erfi(50) 返回 6 * 10 ** 1083。考虑到我对积分的处理，这里有一个精确的表示很好。(2) 我需要做的是计算从 -inf 到 x 的积分，其中 x 在函数调用中从 -150 变化到 50。我不确定该怎么做。 - abcd

不用在意。我最近发表的两个评论实际上是一个关于计算双重积分的单独问题。我打算自己尝试解决这个问题，如果我做不到，就会发布一个新的问题。 - abcd

网页内容由stack overflow 提供, 点击上面的

可以查看英文原文，
原文链接

- ali_m · Accepted Answer

我认为@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))更加数值稳定。

更新

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