Luce規則的python sympy

Question

我想使用Luce's axiom來計算二元決策的概率。相似性函數使用指數如下定義;

sA = b+exp(-|x-xA|) 
sB = b+exp(-|x-xB|) 
pA = 1/(1+(sB/sA)) 

爲了獲得損失，我們需要將pA和1-pA整合到x的各自範圍內。

loss = integrate(pA, (x,0,0.5)) + integrate(1-pA, (x,0.5,1)) 

當使用sympy寫入我得到的損耗當b = 0（0.5075），但以下錯誤當b> 0;

加註PolynomialDivisionFailed（F，G，K） sympy.polys.polyerrors.PolynomialDivisionFailed：將 [-0.426881219248826 * _z + 0.0106631460069606]由[_z時不能降低在一個多項式>除法算法度 - 0.0249791874803424。當係數域中的tect零不可能發生時，這可能發生。計算域是RR（_z）。零檢測>在此係統中保證nt域。這可能表明SymPy中存在錯誤，或者該域是用戶定義的，並且不會正確執行零檢測。

我不確定這個錯誤的含義。

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() 

作爲一個說明工作部分需要幾分鐘的時間來計算，有沒有什麼辦法來加快步伐？

我會感謝任何幫助/建議。

感謝

編輯：編輯代碼

Answer 1

當你真正計算積分一個錯字，你整合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] 

希望這有助於！