# Generalized incomplete gamma function

As I know there is a built in function for incomplete gamma function, incgam(x, a), in R. May I know is there a built in function for the generalized incomplete gamma function? Or how can I modify the incomplete gamma built in function to get the generalized incomplete gamma function?

My generalized incomplete gamma function is $$gamma(a,x;b) = \int_x^\infty t^{a-1}e^{-t-bt^{-1}} dt$$ and the incomplete gamma function is $$gamma(a,x) = \int_x^\infty t^{a-1}e^{-t} dt$$

• 1. Which generalized incomplete gamma do you mean? 2. To the extent that this is a "which function do I call?" question, this is unlikely to be on topic – Glen_b -Reinstate Monica Nov 7 '18 at 0:51
• My generalized incomplete gamma function is $$gamma(a,x;b) = \int_x^\infty t^{a-1}e^{-t-bt^{-1}} dt$$ and the incomplete gamma function is same as provided by Martijn below – Alicia Nov 7 '18 at 1:53
• Thanks; could you edit this crucial information into your question please – Glen_b -Reinstate Monica Nov 7 '18 at 2:14

### For the case $${\Gamma(a,x,y) = \int_x^y t^{a-1}e^{-t} dt}$$

You can compute the generalized incomplete gamma/beta/whatever function as the difference between two upper/lower incomplete functions. $$\begin{array}{rcl}\int_x^y f (t) dt &=& \int_0^y f (t) dt - \int_0^x f (t) dt \\&=& \int_x^1 f (t) dt - \int_y^1 f (t) dt \\&=& \int_x^\infty f (t) dt - \int_y^\infty f (t) dt\end{array}$$

In your case with the incgam function from the pracma package, which is the upper incomplete gamma function $$\int_x^\infty t^{a-1}e^{-t} dt$$ you have the third line and you can program this as

gengam <- function(x, y, a) {
pracma::incgam(x, a) - pracma::incgam(y, a)
}


### For the case $${\Gamma(a,x,b) = \int_x^\infty t^{a-1}e^{-t-bt^{-1}} dt}$$

You can either use

1. a straigtforward integration (e.g. trapezium rule or other type of computation methods)
2. or you use a derived analytical form that might be computed faster.

For this second case there is some closed form solution (in terms of Bessel and error functions) given in Chaudhry, M. Aslam, N. M. Temme, and E. J. M. Veling. "Asymptotics and closed form of a generalized incomplete gamma function." Journal of Computational and Applied Mathematics 67.2 (1996): 371-379.

In the case of half-integer values of $$a$$ ($$\frac{1}{2}, \frac{3}{2}, \frac{5}{2}, etc$$) you can use:

$$\begin{multline} \Gamma(a,x;b) = b^{a/2} \lbrace K_a(2\sqrt{b}) + \pi(-1)^{a-1/2} I_a(2\sqrt{b}) \rbrace \, \text{erfc}(\sqrt{x}+\sqrt{b/x}) \\ + b^{a/2} K_a(2\sqrt{b}) \, \text{erfc}(\sqrt{x}-\sqrt{b/x}) \\+ 2 e^{-x-b/x} \sum_{j=0}^{a-3/2} x^{j+1/2}b^{a/2-1/2j-1/4} \lbrace (-1)^{a+j+1/2} K_{j+1/2}(2\sqrt{b})I_a(2\sqrt{b})\\+I_{j+1/2}(2\sqrt{b})K_a(2\sqrt{b})\rbrace \end{multline}$$

(this is a slightly more simpler special case of an expression where $$a$$ is not restricted to half integer values, see this in Miller, Allen R., and Ira S. Moskowitz. "On certain generalized incomplete gamma functions." Journal of computational and applied mathematics 91.2 (1998): 179-190 )

Comparison of the two methods in R

genincgam <- function(x, a, b) {
if (mod(a+0.5,1) != 0) {
stop("the variable 'a' must be a half integer")
}

# calculating pieces
# storing K1 and I1 which are used several times
K1 <- besselK(2*sqrt(b),a)
I1 <- besselI(2*sqrt(b),a)
# summation term
j <- c(0:(max(0,a-1.5)))
S <- sum( x^(j+0.5) * b^(a/2-0.5*j-0.25) * ( (-1)^(a+j+0.5) * besselK(2*sqrt(b),j+0.5)*I1 + besselI(2*sqrt(b),j+0.5)*K1))

b^(a/2)*(K1 + pi *(-1)^(a-0.5)*I1) * erfc(sqrt(x)+sqrt(b/x)) +
b^(a/2)*K1*erfc(sqrt(x)-sqrt(b/x)) +
2*exp(-x-b/x) * S
}

f <- function(t, a, b) {
t^(a-1)*exp(-t-b/t)
}

> # method 2
> genincgam(2, a = 0.5, b = 1)
[1] 0.05573681
> # method 1
> integrate(g <- function(t) f(t,a = 0.5, b = 1),  2, Inf, rel.tol=10^{-8})
0.05573681 with absolute error < 2.4e-09
> integrate(g <- function(t) f(t,a = 0.5, b = 1),  2, Inf, rel.tol=10^{-5})
0.05573681 with absolute error < 7.5e-07
> integrate(g <- function(t) f(t,a = 0.5, b = 1),  2, Inf, rel.tol=10^{-2})
0.05573628 with absolute error < 0.001
> #compare computation time
> system.time( replicate(10000,genincgam(2, a = 0.5, b = 1) ) )
user  system elapsed
0.312   0.000   0.312
> system.time( replicate(10000,integrate(g <- function(t) f(t,a = 0.5, b = 1),  2, Inf, rel.tol=10^{-8}) ) )
user  system elapsed
0.788   0.000   0.787
> system.time( replicate(10000,integrate(g <- function(t) f(t,a = 0.5, b = 1),  2, Inf, rel.tol=10^{-5}) ) )
user  system elapsed
0.628   0.000   0.628
> system.time( replicate(10000,integrate(g <- function(t) f(t,a = 0.5, b = 1),  2, Inf, rel.tol=10^{-2}) ) )
user  system elapsed
0.3     0.0     0.3


It will depend somewhat on the values of the variables and the desired precision but in general the second case will be faster.

• Hi Martijn, thank you for your reply. May I know how it going to put in R by using the incomplete gamma function, incgam(x, a)? – Alicia Nov 6 '18 at 9:01
• Actually I have a cdf which having the generalized incomplete gamma function such as in the form of (v/p, ax^p; ab). (Sorry, I don't know how to use the mathematics code.) – Alicia Nov 6 '18 at 10:18
• but my generalized incomplete gamma function is $$\int_x^\infty t^{a-1}e^{-t-bt^{-1}} dt$$. Is it still using your suggest way to find it? – Alicia Nov 7 '18 at 1:42
• $$gamma(a,x;b) = \int_x^\infty t^{a-1}e^{-t-bt^{-1}} dt$$ – Alicia Nov 7 '18 at 1:49
• Ah I see. I should have anticipated that since generalisation is an ambiguous thing. – Sextus Empiricus Nov 7 '18 at 7:33
The text passage is a discussion about the generalized incomplete gamma function, including its definition and methods to compute it in R. It also includes a comparison of two methods for computing a specific case of the function.
no_question

the passage does not contain a direct mathematical question that requires an answer. instead, it provides information and explanations about the generalized incomplete gamma function and how to compute it in r.
Given a text passage, determine:
    1. If it contains a mathematical question (including word problems or calculations)
    2. If an answer to that mathematical question is provided in the text

First, describe the content structure you found in 1-2 sentences.

Then provide your analysis in the following format:
TYPE: <one of the following values>
- NO_QUESTION
- SINGLE_QUESTION
- SINGLE_QUESTION_WITH_ANSWER
- MULTIPLE_QUESTIONS
- MULTIPLE_QUESTIONS_WITH_ANSWERS
system
# Generalized incomplete gamma function

As I know there is a built in function for incomplete gamma function, incgam(x, a), in R. May I know is there a built in function for the generalized incomplete gamma function? Or how can I modify the incomplete gamma built in function to get the generalized incomplete gamma function?

My generalized incomplete gamma function is $$gamma(a,x;b) = \int_x^\infty t^{a-1}e^{-t-bt^{-1}} dt$$ and the incomplete gamma function is $$gamma(a,x) = \int_x^\infty t^{a-1}e^{-t} dt$$

• 1. Which generalized incomplete gamma do you mean? 2. To the extent that this is a "which function do I call?" question, this is unlikely to be on topic – Glen_b -Reinstate Monica Nov 7 '18 at 0:51
• My generalized incomplete gamma function is $$gamma(a,x;b) = \int_x^\infty t^{a-1}e^{-t-bt^{-1}} dt$$ and the incomplete gamma function is same as provided by Martijn below – Alicia Nov 7 '18 at 1:53
• Thanks; could you edit this crucial information into your question please – Glen_b -Reinstate Monica Nov 7 '18 at 2:14

### For the case $${\Gamma(a,x,y) = \int_x^y t^{a-1}e^{-t} dt}$$

You can compute the generalized incomplete gamma/beta/whatever function as the difference between two upper/lower incomplete functions. $$\begin{array}{rcl}\int_x^y f (t) dt &=& \int_0^y f (t) dt - \int_0^x f (t) dt \\&=& \int_x^1 f (t) dt - \int_y^1 f (t) dt \\&=& \int_x^\infty f (t) dt - \int_y^\infty f (t) dt\end{array}$$

In your case with the incgam function from the pracma package, which is the upper incomplete gamma function $$\int_x^\infty t^{a-1}e^{-t} dt$$ you have the third line and you can program this as

gengam <- function(x, y, a) {
pracma::incgam(x, a) - pracma::incgam(y, a)
}


### For the case $${\Gamma(a,x,b) = \int_x^\infty t^{a-1}e^{-t-bt^{-1}} dt}$$

You can either use

1. a straigtforward integration (e.g. trapezium rule or other type of computation methods)
2. or you use a derived analytical form that might be computed faster.

For this second case there is some closed form solution (in terms of Bessel and error functions) given in Chaudhry, M. Aslam, N. M. Temme, and E. J. M. Veling. "Asymptotics and closed form of a generalized incomplete gamma function." Journal of Computational and Applied Mathematics 67.2 (1996): 371-379.

In the case of half-integer values of $$a$$ ($$\frac{1}{2}, \frac{3}{2}, \frac{5}{2}, etc$$) you can use:

$$\begin{multline} \Gamma(a,x;b) = b^{a/2} \lbrace K_a(2\sqrt{b}) + \pi(-1)^{a-1/2} I_a(2\sqrt{b}) \rbrace \, \text{erfc}(\sqrt{x}+\sqrt{b/x}) \\ + b^{a/2} K_a(2\sqrt{b}) \, \text{erfc}(\sqrt{x}-\sqrt{b/x}) \\+ 2 e^{-x-b/x} \sum_{j=0}^{a-3/2} x^{j+1/2}b^{a/2-1/2j-1/4} \lbrace (-1)^{a+j+1/2} K_{j+1/2}(2\sqrt{b})I_a(2\sqrt{b})\\+I_{j+1/2}(2\sqrt{b})K_a(2\sqrt{b})\rbrace \end{multline}$$

(this is a slightly more simpler special case of an expression where $$a$$ is not restricted to half integer values, see this in Miller, Allen R., and Ira S. Moskowitz. "On certain generalized incomplete gamma functions." Journal of computational and applied mathematics 91.2 (1998): 179-190 )

Comparison of the two methods in R

genincgam <- function(x, a, b) {
if (mod(a+0.5,1) != 0) {
stop("the variable 'a' must be a half integer")
}

# calculating pieces
# storing K1 and I1 which are used several times
K1 <- besselK(2*sqrt(b),a)
I1 <- besselI(2*sqrt(b),a)
# summation term
j <- c(0:(max(0,a-1.5)))
S <- sum( x^(j+0.5) * b^(a/2-0.5*j-0.25) * ( (-1)^(a+j+0.5) * besselK(2*sqrt(b),j+0.5)*I1 + besselI(2*sqrt(b),j+0.5)*K1))

b^(a/2)*(K1 + pi *(-1)^(a-0.5)*I1) * erfc(sqrt(x)+sqrt(b/x)) +
b^(a/2)*K1*erfc(sqrt(x)-sqrt(b/x)) +
2*exp(-x-b/x) * S
}

f <- function(t, a, b) {
t^(a-1)*exp(-t-b/t)
}

> # method 2
> genincgam(2, a = 0.5, b = 1)
[1] 0.05573681
> # method 1
> integrate(g <- function(t) f(t,a = 0.5, b = 1),  2, Inf, rel.tol=10^{-8})
0.05573681 with absolute error < 2.4e-09
> integrate(g <- function(t) f(t,a = 0.5, b = 1),  2, Inf, rel.tol=10^{-5})
0.05573681 with absolute error < 7.5e-07
> integrate(g <- function(t) f(t,a = 0.5, b = 1),  2, Inf, rel.tol=10^{-2})
0.05573628 with absolute error < 0.001
> #compare computation time
> system.time( replicate(10000,genincgam(2, a = 0.5, b = 1) ) )
user  system elapsed
0.312   0.000   0.312
> system.time( replicate(10000,integrate(g <- function(t) f(t,a = 0.5, b = 1),  2, Inf, rel.tol=10^{-8}) ) )
user  system elapsed
0.788   0.000   0.787
> system.time( replicate(10000,integrate(g <- function(t) f(t,a = 0.5, b = 1),  2, Inf, rel.tol=10^{-5}) ) )
user  system elapsed
0.628   0.000   0.628
> system.time( replicate(10000,integrate(g <- function(t) f(t,a = 0.5, b = 1),  2, Inf, rel.tol=10^{-2}) ) )
user  system elapsed
0.3     0.0     0.3


It will depend somewhat on the values of the variables and the desired precision but in general the second case will be faster.

• Hi Martijn, thank you for your reply. May I know how it going to put in R by using the incomplete gamma function, incgam(x, a)? – Alicia Nov 6 '18 at 9:01
• Actually I have a cdf which having the generalized incomplete gamma function such as in the form of (v/p, ax^p; ab). (Sorry, I don't know how to use the mathematics code.) – Alicia Nov 6 '18 at 10:18
• but my generalized incomplete gamma function is $$\int_x^\infty t^{a-1}e^{-t-bt^{-1}} dt$$. Is it still using your suggest way to find it? – Alicia Nov 7 '18 at 1:42
• $$gamma(a,x;b) = \int_x^\infty t^{a-1}e^{-t-bt^{-1}} dt$$ – Alicia Nov 7 '18 at 1:49
• Ah I see. I should have anticipated that since generalisation is an ambiguous thing. – Sextus Empiricus Nov 7 '18 at 7:33
user
The text passage is a discussion about the generalized incomplete gamma function, including its definition and methods to compute it in R. It also includes a comparison of two methods for computing a specific case of the function.

TYPE: NO_QUESTION

The passage does not contain a direct mathematical question that requires an answer. Instead, it provides information and explanations about the generalized incomplete gamma function and how to compute it in R.