I want to calculate probability for a binary decision using Luce's axiom. The similarity function is defined using exponential as follows;
sA = b+exp(-|x-xA|)
sB = b+exp(-|x-xB|)
pA = 1/(1+(sB/sA))
To get the loss we need to integrate both pA and 1-pA over x in their respective ranges.
loss = integrate(pA, (x,0,0.5)) + integrate(1-pA, (x,0.5,1))
When written using sympy I get the loss when b=0 (0.5075) but the following error when 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.
I am not sure what this error means.
The python code is (error does not depend on specific xA and 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()
As a note the working part takes a few minutes to compute, is there any way to speed it up?
When you actually evaluate the integral, you integrate pA in the second integral.
In the description you say it should be 1 - pA, so I'll assume that's what you want.
The fact that the integral doesn't evaluate appears to be a bug in SymPy. Here is a modification that works on my machine.
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)))
Unfortunately this is still horribly slow. Both the fact that this works and the fact that it takes so long could be idiosyncrasies of my installation since I install the development version of sympy regularly.
I would really suggest using some form of numerical integration unless you specifically need a symbolic expression. Given the same initialization and imports above (but not the integral) this can be done like this:
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 also has some integration routines. The following would be an alternative to the above two lines.
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]
Hope this helps!