python sympy for Luce's rule

Question

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?

I will appreciate any help/suggestions.

Thanks

EDIT: edited a typo in the code

Answer 1

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!