在 exp 的标量除法和溢出中遇到 numpy 无效值

Numpy invalid value encountered in scalar divide and overflow in exp

提问人:jasmine 提问时间:11/8/2023 最后编辑:jasmine 更新时间:11/8/2023 访问量:121

问:

这是我试图计算的:

    def attraction(b1, b2, b3):
        a1=10 * b1 + 20 * b2 + 6 * b3
        a2=12 * b1 + 18 * b2 + 10 * b3
        a3=0 * b1 + 10 * b2 + 30 * b3
        return a1, a2, a3

    def prob(a1, a2, a3, phi):
        de = np.exp(phi * (a1 - a2)) + 1 + np.exp(phi*(a3-a2))
        p1 = np.exp(phi * (a1 - a2)) / de
        p2 = 1 / de
        p3 = np.exp(phi * (a3 - a2)) / de
        return p1, p2, p3

b1, b2, b3是概率,所以在和之间,在和之间,在和之间。 我正在做优化,是一个要估计的参数,它没有边界,但它从来没有给我一个大于 1 或小于 0 的值。a1620a21218a3030phi

我有时会得到

RuntimeWarning: overflow encountered in exp p3 = np.exp(phi * (a3 - a2)) / de

RuntimeWarning: invalid value encountered in scalar divide p3 = np.exp(phi * (a3 - a2)) / de

RuntimeWarning: overflow encountered in exp
  de = np.exp(phi * (a1 - a2)) + 1 + np.exp(phi*(a3-a2))

我想这不是因为是零,而是因为太大了。 我尝试设置,但出现错误。经过一番搜索,我想在 M1 mac 中不可用。dedefloat128TypeError: data type 'float128' not understoodfloat128

我的问题是如何在不放弃精度的情况下解决问题。

python numpy apple-m1 整数溢出

评论

1赞 Sembei Norimaki 11/8/2023
请将代码发布到您定义参数值的位置。
0赞 lastchance 11/8/2023
您确定 a1、a2、a3 和 phi 的这些限制吗?您将得到的最大术语是 exp(30),这不应该溢出。
0赞 jasmine 11/8/2023
@lastchance确切地说,这就是我感到困惑的原因。
0赞 lastchance 11/8/2023
喇嘛达从何而来?它会停止 p1、p2 和 p3 加起来为 1(并且与您的错误消息不对应)。
0赞 jasmine 11/8/2023
@lastchance是的,对不起,我复制了错误的版本。它应该是 phi,已更正。

答:

0赞 Kunal Shah 11/8/2023 #1

您可以修改函数以使用对数总和而不是指数。此链接提供了一些很好的例子,并详细介绍了数学基础。

来自文章:

>>> x = np.array([1000, 1000, 1000])
>>> np.exp(x)
array([inf, inf, inf])

即使指数为 1000,该值也会溢出。然后,他们添加一个函数:

def logsumexp(x):
    c = x.max()
    return c + np.log(np.sum(np.exp(x - c)))

这使您可以间接计算指数

>>> logsumexp(x)
1001.0986122886682
>>> np.exp(x - logsumexp(x))
array([0.33333333, 0.33333333, 0.33333333])
0赞 ti7 11/8/2023 #2

我无法使用给定范围内的值重现这一点(也许是非常大或 1
但是,您可以使用 SymPy 解决这个问题,因为它可以以非常高的精度表示值,并且既可以与 NumPy 的类型限制一起使用,也可以逐步消除 NumPy 的类型限制!(**咳嗽,执行速度受到巨大惩罚)
lamdanan

设置

>>> from sympy import symbols, exp
>>> a1, a2, a3, phi = symbols("a1 a2 a3 phi")
>>> de = exp(phi * (a1 - a2)) + 1 + exp(phi*(a3-a2))
>>> p1 = exp(phi * (a1 - a2)) / de
>>> p1
exp(phi*(a1 - a2))/(exp(phi*(a1 - a2)) + exp(phi*(-a2 + a3)) + 1)

Lambdify 创建一个适用于 NumPy 类型的函数

>>> expr = lambdify((a1, a2, a3, phi), p1)

试一试

>>> a1_ = np.linspace( 6, 20, 20)
>>> a2_ = np.linspace(12, 18, 20)
>>> a3_ = np.linspace( 0, 30, 20)
>>> phx = np.linspace( 0,  1, 20)
>>> expr(a1=a1_, a2=a2_, a3=a3_, phi=phx)
array([3.33333333e-01, 3.22212128e-01, 2.97978285e-01, 2.70970702e-01,
       2.46581616e-01, 2.26688815e-01, 2.11128138e-01, 1.98368944e-01,
       1.85454387e-01, 1.67730136e-01, 1.40043425e-01, 1.01808120e-01,
       6.17809832e-02, 3.12962524e-02, 1.36403988e-02, 5.27442099e-03,
       1.84280490e-03, 5.86558449e-04, 1.70613612e-04, 4.53975898e-05])

原始功能(请注意,这不会重现,但我没有 M1 拱门系统,也没有lamda)

>>> lamda = 1
>>> prob(a1_, a2_, a3_, phx)[0]
array([3.33333333e-01, 3.22212128e-01, 2.97978285e-01, 2.70970702e-01,
       2.46581616e-01, 2.26688815e-01, 2.11128138e-01, 1.98368944e-01,
       1.85454387e-01, 1.67730136e-01, 1.40043425e-01, 1.01808120e-01,
       6.17809832e-02, 3.12962524e-02, 1.36403988e-02, 5.27442099e-03,
       1.84280490e-03, 5.86558449e-04, 1.70613612e-04, 4.53975898e-05])

SymPy 也可以简化表达式,但不会太多

>>> p1.simplify()
exp(a1*phi)/(exp(a1*phi) + exp(a2*phi) + exp(a3*phi))

一些尺寸的演示

>>> np.exp(1000)
inf
>>> sympy.exp(1000)
exp(1000)
>>> sympy.exp(1000).n()
1.97007111401705e+434
>>> exp(10**6) / pi
exp(1000000)/pi
>>> (exp(10**6) / pi).n()
9.65502447726994e+434293
>>> sympy.exp(1000).n() / np.arange(1,5)  # note 2
array([1.97007111401705e+434, 9.85035557008524e+433,
       6.56690371339016e+433, 4.92517778504262e+433], dtype=object)

脚注

  1. 评论有关导致/包含已删除的问题denan
  2. 请注意,这里的数组内容不是一个正确的数字 NumPy 数组,而是包含对象,这些对象是 Python 类型,使用效率低几个数量级