我想使用 Luce's axiom 计算二元决策的概率。相似度函数使用指数定义如下;
sA = b+exp(-|x-xA|)
sB = b+exp(-|x-xB|)
pA = 1/(1+(sB/sA))
为了获得损耗,我们需要在各自的范围内对 x 上的 pA 和 1-pA 进行积分。
loss = integrate(pA, (x,0,0.5)) + integrate(1-pA, (x,0.5,1))
当使用 sympy 编写时,当 b=0 (0.5075) 时我得到损失,但当 b>0 时出现以下错误;
raise PolynomialDivisionFailed(f, g, K) sympy.polys.polyerrors.PolynomialDivisionFailed: couldn't reduce degree in a polynomial >division algorithm when dividing [-0.426881219248826*_z + 0.0106631460069606] by [_z - 0.0249791874803424]. This can >happen when it's not possible to de tect zero in the coefficient domain. The domain of computation is RR(_z). Zero detection >is guaranteed in this coefficie nt domain. This may indicate a bug in SymPy or the domain is user defined and doesn't >implement zero detection properly.
我不确定这个错误意味着什么。
Python代码为(错误不依赖于具体的xA和xB);
from sympy import *
var('x')
xA = 0.8
xB = 0.9
#this works
b = 0
sA = b+exp(-abs(x-xA))
sB = b+exp(-abs(x-xB))
pA = 1/(1+(sB/sA))
print pA
loss = integrate(pA, (x,0,0.5)) + integrate(1-pA, (x,0.5,1))
print loss.evalf()
#this doesn't
b = 1
sA = b+exp(-abs(x-xA))
sB = b+exp(-abs(x-xB))
pA = 1/(1+(sB/sA))
print pA
loss = integrate(pA, (x,0,0.5)) + integrate(1-pA, (x,0.5,1)) #fails here
print loss.evalf()
请注意,工作部分需要几分钟才能计算,有什么方法可以加快速度吗?
我将不胜感激任何帮助/建议。
谢谢
编辑:编辑了代码中的拼写错误
最佳答案
当您实际计算积分时,您将 pA
积分到第二个积分中。
在描述中,您说它应该是 1 - pA
,所以我假设这就是您想要的。
积分不计算的事实似乎是 SymPy 中的一个错误。 这是在我的机器上运行的修改。
import sympy as sy
x = sy.symbols('x')
b = 1
sA = b + sy.exp(- sy.Abs(x - xA))
sB = b + sy.exp(- sy.Abs(x - xB))
pA = 1 / (1 + (sB / sA))
sy.N(sy.Integral(pA, (x, 0, 0.5)) + sy.Integral(1 - pA, (x, 0.5, 1)))
不幸的是,这仍然非常慢。 由于我定期安装 sympy 的开发版本,因此它有效以及花费很长时间的事实可能是我的安装特性。
我真的建议使用某种形式的数值积分,除非您特别需要符号表达式。 给定与上面相同的初始化和导入(但不是积分),可以这样完成:
from sympy.mpmath import quad
# Make the expression into a callable function.
pA_func = sy.lambdify([x], pA)
quad(pA_func, [0, .5]) + quad(lambda x: 1 - pA_func(x), [.5, 1])
SciPy 也有一些集成例程。 以下是上述两行的替代。
from scipy.integrate import quad
# Make the expression into a callable function.
pA_func = sy.lambdify([x], pA)
quad(pA_func, 0, .5)[0] + quad(lambda x: 1 - pA_func(x), .5, 1)[0]
希望这有帮助!
关于卢斯规则的 python sympy,我们在Stack Overflow上找到一个类似的问题: https://stackoverflow.com/questions/21353081/