4
$\begingroup$

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$$

$\endgroup$
  • $\begingroup$ 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 $\endgroup$ – Glen_b -Reinstate Monica Nov 7 '18 at 0:51
  • $\begingroup$ 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 $\endgroup$ – Alicia Nov 7 '18 at 1:53
  • $\begingroup$ Thanks; could you edit this crucial information into your question please $\endgroup$ – Glen_b -Reinstate Monica Nov 7 '18 at 2:14
4
$\begingroup$

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.

$\endgroup$
  • $\begingroup$ 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)? $\endgroup$ – Alicia Nov 6 '18 at 9:01
  • $\begingroup$ 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.) $\endgroup$ – Alicia Nov 6 '18 at 10:18
  • $\begingroup$ 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? $\endgroup$ – Alicia Nov 7 '18 at 1:42
  • $\begingroup$ $$gamma(a,x;b) = \int_x^\infty t^{a-1}e^{-t-bt^{-1}} dt$$ $\endgroup$ – Alicia Nov 7 '18 at 1:49
  • $\begingroup$ Ah I see. I should have anticipated that since generalisation is an ambiguous thing. $\endgroup$ – Sextus Empiricus Nov 7 '18 at 7:33

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.