2

I have used SageMath to symbolically integrate an expression. The result contains a gamma function with two input parameters:

gamma(-1, 2*((x - xp)^2 + (y - yp)^2)/s^2)

Apparently this is called an incomplete gamma function. Now I want to use this expression in a Python code. I have tracked down the incomplete gamma function to scipy.special.gammainc. Unfortunately, this function does not allow negative input parameters and I have to use -1 as the first input parameter. How can I work around this?

SiHa
  • 7,830
  • 13
  • 34
  • 43
Numaerius
  • 383
  • 3
  • 10
  • the gamma function [isn't defined for negative integers](https://en.wikipedia.org/wiki/Gamma_function). Why do you need **-1** as a parameter? – DrBwts Jan 15 '20 at 16:17
  • That is the expression that SageMath spits out. SageMath can evaluate `gamma(-1, someValue)` and approximate this numerically. – Numaerius Jan 15 '20 at 16:19
  • I've rolled back your edit. It is perfectly normal to answer your own question, and you have done so in the appropriate place. There is no need to draw attention to the answer in the question itself. Note that you are also free to accept your own answer (but probably have to wait for a grace period first) – SiHa Jan 16 '20 at 13:16
  • Thanks @SiHa. I indeed have to wait for a day before I can accept my own answer. – Numaerius Jan 16 '20 at 13:18

2 Answers2

3

The lower incomplete gamma function can be defined in terms of an improper integral according to Wikipedia. This integral can be related to the generalised form of the exponential integral. Both pages give this relation between the two:

E_n(x) = x^(n-1)*gamma(1-n, x)

So for the case in OP we would have n=2, which corresponds to -1 as a first input parameter for the gamma function. I have numerically verified in SageMath that the above relation holds up. The corresponding functions in SageMath are:

1. gamma(n, x) == gamma_inc(n, x)
2. E_n(x) == exp_integral_e(n, x)

Which according to the Wikipedia relation gives us (aside from round-off errors):

exp_integral_e(n, x) == x^(n-1)*gamma_inc(1-n, x)

The corresponding Python functions are:

1. gamma(n, x) == scipy.special.gammainc(n, x)
2. E_n(x) == scipy.special.expn(n, x)

Which according to the Wikipedia relation gives us (aside from round-off errors):

expn(n, x) == x**(n-1)*gammainc(1-n, x)

There is one small caveat. The gamma_inc function from SageMath accepts a negative first input parameter, whereas the gammainc function from scipy.special does not. However, the expn function from scipy.special does not have this limitation as it can be evaluated for n>=2 corresponding to a negative first input parameter for gamma_inc.

So the answer to the OP is to use the Wikipedia relation to replace the lower incomplete gamma function with the generalised exponential integral and to use scipy.special.expn for evaluation in Python.

Numaerius
  • 383
  • 3
  • 10
-2

There is a reason you cannot input a negative number, as the factorials of negative integers cannot be computed, since for n = 0, the recurrence relation is: 

(n-1)! = n!/n

This would result in division by zero.
Perhaps you should rephrase your question into what type of goal you are trying to accomplish, then cross post https://math.stackexchange.com/

tbrk
  • 156
  • 8
  • 1
    I am not trying to argue with the math and I am also not looking for a mathematical proof. I am merely trying to locate a Python function which can do the same as the `gamma` function from SageMath. Whether it should be possible or not, the SageMath solution with `gamma(-1, someValue)` works as intended. – Numaerius Jan 15 '20 at 17:28