About incomplete gamma function

Hi,

I require the incomplete gamma function for my MLE calculations. I realize that the implementations in Mathematica and matlab are different.

Basically I want to calculate value of this integral,

int(t^(a-1)*exp(-t), t, x, inf), where a and x are inputs to the incomplete gamma function.

The problem here is not that I want the upper tail, but that matlab implementation which is,

gammainc(x,a), is different from what I want i.e. gammainc(a,x), the parameters in this case are reversed.

But matlab does not allow negative value for a, which is why a simple swapping of input parameters is not working in my case.

If I try to integrate using matlab's int function, the speed of execution is very slow (though I get the desired answer). I was wondering if someone could suggest a faster way for me to do this. My guess is that matlab's gammainc uses mex files to speed things up, but if someone could suggest a more convenient way, I would really appreciate it.

Thanks,

Yogesh

Answer 1

[Mike Hosea]

I'm confused by the swapping parameters aspect, perhaps because I have not used Mathematica for nearly two decades. Could we forget about gammainc functions in MATLAB and Mathematica for a moment and just tell me the integral you want to evaluate and the ranges of the two parameters that you want it to work for? For example, "I want to evaluate the integral of t^(a-1)*exp(-t) from x to infinity where a > 0 and x >= 1" or whatever it has to be. -- Mike

[YOGESH VIRKAR]

Hi Mike,

I want to integrate,

int(t^(a-1)*exp(-t), t, x, inf), where x,a are real. Actually, this integration is taking place in a function which I am trying to minimize using fminsearch. The integration is a part of MLE calculation.

Now the problem I am facing here is of speed, since the parameter x has about a 100 values, which means I need to evaluate the integral 100 times. Also fminsearch would be repeating these calculations till it finds the minimum, which makes the process really slow.

Thanks,

Yogesh

[Walter Roberson]

That integral is the "upper tail" of the incomplete gamma function. It is, however, not defined if a is 0 or a negative integer.

Earlier today I explored (in Maple) a number of ways of expressing it, but I will need to go back and revisit them. The most compact way I found involved hypergeom() with [a+2] as its second parameter; unfortunately when a < -2 that leaves a negative value for the second parameter and Maple's hypergeom() complains about that. I may need to rewrite it in terms of pochhammer.

Maple had no difficulty evaluating GAMMA(a,x) for negative values of "a"; Maple's GAMMA function is equivalent to the "upper tail", exactly the value you were looking to integrate (except perhaps to within a scale factor of 1/(a+1) )

[YOGESH VIRKAR]

Well, even Mathematica has a Gamma[a,x] function which is exactly the integral I am looking for and it does allow negative values for a.

One thing I did not understand is that if for the upper tail of the gamma function 'a' should be non-negative, why does directly evaluating the integral (for upper tail of gamma function) gives me the desired result? That is why, I was confused, I thought on doing such an integration, using int function in matlab I might get Inf or NaN, but the answer I get agrees with the one I get from Mathematica's Gamma[a,x].

The only problem with int function in matlab is that of speed, which makes the approach not feasible for my application.

[Mike Hosea]

If you will allow the constraint x > 0, you can use QUADGK (with whatever tolerances you like). You have to put it in a loop for vector inputs:

function y = foo(x,a)

y = zeros(size(x));

for k = 1:numel(y)

y(k) = quadgk(@(t)t.^(a(k)-1).*exp(-t),x(k),inf,'abstol',1e-12,'reltol',1e-12);

end

[YOGESH VIRKAR]

Hi Mike,

My x's are always positive....thanks a ton. That works! I had this big deadline in two days, and now I am good!

Thanks,

Yogesh