# Numerical Double integration with endpoint singularity in scipy Python gives incorrect answer

I am trying to integrate the following function in Python,

$\int_{0}^{\infty}\int_{0}^{\infty} \dfrac{e^{-x-y}}{B(x,y)}dx dy$,

where $B(x,y)$ is the beta function - $B(x,y) = \int_0^{1}a^{x-1}(1-a)^{y-1}da$.

As the above function was not directly integrable in scipy.integrate, I did a change of variable, $x\rightarrow e^{-x},y\rightarrow e^{-y}$ to have

$\int_{0}^{1}\int_{0}^{1} \dfrac{1}{B(-\log (x),-\log (y))}dx dy$.

scipy.integrate.dblquad(lambda y,x: 1/ scipy.special.beta(-math.log(x),-math.log(y) ), 0, 1, lambda x: 0, lambda x: 1)

returns $6.81\ldots$ with error in order of $10^{-6}$. But integrating the above functions (with and without change of variable) in Mathematica (using NIntegrate) returns $2.0505\ldots$. I do not understand this huge difference in answer when python returns an answer with $10^{-6}$ accuracy. Whats happening? Is there a better way to integrate it in python? I plan to perform the integration many number of times as part of an algorithm written in python.

$\frac{1}{B(-\log (x),-\log (y))}$ looks like this Edit

Can someone with Windows 10 and Anaconda reproduce the results I get?

I get the incorrect answer using Windows 10, Python 3.5.2 |Anaconda 4.3.0 (64-bit), scipy version 0.18.1. I am not sure about the BLAS linkage.

I get the correct answer in Linux using Python 2.7.12, scipy version 0.18.1

• I get (2.0505570600455982, 4.867808183917077e-05) when I run the python code (scipy version: 0.18.1). Feb 12 '17 at 20:19
• Wow! Mine is also 0.18.1, using Python 3.5 in Spyder IDE. I don't get it..Though I get this warning IntegrationWarning: The integral is probably divergent, or slowly convergent. Feb 12 '17 at 20:23
• @ Stelios - I tried it in an online python here link and I get (2.0505570600455982, 4.867808183917077e-05) . But in Ipython console in Spyder in my local machine I get (6.8194728123941895, 2.7699807105818763e-06) . I tried it in a friends laptop with Spyder and I get 6.81. as well. I am confused.. Maybe Spyder has something wrong.. Feb 12 '17 at 20:52
• Very weird... I use Python 3.6 and Jupyter (IPython), BTW. Try posting this on stackoverflow SE to see if this issue can be reproduced by other users. Feb 12 '17 at 21:17
• Yes, My issue was also resolved by reinstalling numpy and scipy from lfd.uci.edu/~gohlke/pythonlibs/#scipy .Thanks Krill for pointing it out! Any idea how to send a bug report to anaconda? Feb 14 '17 at 22:51