提问人:jasmine 提问时间:11/8/2023 最后编辑:jasmine 更新时间:11/8/2023 访问量:121
在 exp 的标量除法和溢出中遇到 numpy 无效值
Numpy invalid value encountered in scalar divide and overflow in exp
问:
这是我试图计算的:
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 的值。a1
6
20
a2
12
18
a3
0
30
phi
我有时会得到
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 中不可用。de
de
float128
TypeError: data type 'float128' not understood
float128
我的问题是如何在不放弃精度的情况下解决问题。
答:
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 的类型限制!(**咳嗽,执行速度受到巨大惩罚)lamda
nan
设置
>>> 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)
脚注
- 评论有关导致/包含已删除的问题
de
nan
- 请注意,这里的数组内容不是一个正确的数字 NumPy 数组,而是包含对象,这些对象是 Python 类型,使用效率低几个数量级
评论