I was playing with the Gamma Distribution Function in python. As typically, I used a combination of
bumpy numpy arrays and the scipy.special functions for Gamma and incomplete gamma functions. At one point I realized that the Gamma distribution function (bottom panel on the figure below) was not reaching unity for large x values and for some parameters (k=0.63, theta=0.05). In the given case, the deviation was about 0.3 in probability space, which is significant.
The solution I came up with is to use the module mpmath. When I do, everything looks fine.
I googled only briefly, but did not find any description of this error in scipy.special. Also I am not sure, if using mpmath is the best/fastest/most elegant solution. However, at least it seems to work.