[R] Incomplete Gamma function

[(Ted Harding)]

On 31-Aug-07 13:06:42, Prof Brian Ripley wrote:
> On Fri, 31 Aug 2007, Robin Hankin wrote:
> 
>> Hi Kris
>> lgamma() gives the log of the gamma function.
> 
> Yes, but he used Igamma.  According to ?pgamma,
> 
>   'pgamma' is closely related to the incomplete gamma function.
>   As defined by Abramowitz and Stegun 6.5.1
> 
>      P(a,x) = 1/Gamma(a) integral_0^x t^(a-1) exp(-t) dt
> 
>   P(a, x) is 'pgamma(x, a)'.  Other authors (for example Karl
>   Pearson in his 1922 tables) omit the normalizing factor,
>   defining the incomplete gamma function as
>   'pgamma(x, a) * gamma(a)'.
> 
> and that seems to be what Igamma is following. GSL on the other
> hand has the other tail, so
> 
>> a <- 9
>> x <- 11.1
>> pgamma(x, a, lower=FALSE)*gamma(a)
> [1] 9000.501
> 
>> You need gamma_inc() of the gsl package, a wrapper for the
>> GSL library:
>>
>> > gamma_inc(9,11.1)
>> [1] 9000.501
>> >
> 
> As the above shows, you don't *need* this, but you do need the GSL 
> documentation to find out what R package gsl does.  Why it differs from
> the usual references is something for you to explain.  Wikipedia
> http://en.wikipedia.org/wiki/Incomplete_gamma_function
> distinguishes them, as does MathWorld.
> 
> I suggest you add a clarification to the gsl package as to what the 
> 'incomplete gamma function' means there.

We have been here before! -- though in connection with the
Beta function in the first instance. See:

See the thread starting on 13 Dec 2005 at
  http://finzi.psych.upenn.edu/R/Rhelp02a/archive/66670.html

In particular I'll repeat my views on the distinction of
terminology between "Incomplete Beta/Gamma Function" and
"Beta/Gamma Distribution" (where "Function" refers to the
incomplete *integral* and "Distribution" to the same divided
by the complete integral i.e. by the "Beta/Gamma Function"
which, in my view, should be defined as the complete integral).

  My reasons for preferring the terminology "Incomplete
  ... Function" for the incomplete integral *not* divided
  by the normalising constant (for both Beta and Gamma),
  and using "Distribution" for the incomplete integral
  divided by the constant (i.e. Pearson's "Ratio"), are
  several, but in summary:

  1. The Beta and Gamma functions (not normalised) are
     fundamental mathematical functions in their own right;
     likewise their incomplete versions.

  2. When needed in probability applications, then of course
     they need to be normalised; but then why not simply call
     them "distributions"?

  3. (1) and (2) encapsulate in the terminology an essential
     distinction, and using (2) instead of (1) could lead to
     interesting inferences (e.g. that the complete Beta
     function is identically 1).

     I.e. the Beta function should not change its definition
     as x passes from 1 - epsilon to 1. And similarly for
     the Gamma.

  Granted there is non-uniformity of usage; but this does
  lead to confusion, which could be avoided by simply sticking
  to the distinction between "Incomplete ... Function" and "...
   Distribution". 

For more detail, see
  http://finzi.psych.upenn.edu/R/Rhelp02a/archive/66717.html
(where, in particular, it is pointed out that both Karl Pearson
and Abramowitz and Stegun are inconsistent, within the same
publication, in their terminology, using "... Function" in one
place to mean the integral, in another to mean the probability
distribution. So it is unwise to appeal to either as the
definitive reference, since the outcome will depend on where
in the book you look it up).

It looks as though the documentation for Igamma (ZipfR package)
at

  http://finzi.psych.upenn.edu/R/library/zipfR/html/beta_gamma.html

is admirably explicit as to how this (and related functions) are
defined, so in this case there is no ambiguity.

In the documenation for the Gamma functions in the gsl package,
it is simply stated

  All functions [including gamma_inc()] as documented in the
  GSL reference manual section 7.19.

There is no function named "gamma_inc" in the GSL reference
manual. See:

http://www.gnu.org/software/gsl/manual/html_node/Function-Index.html

All functions are named like "gsl_sf_gamma_inc", so
presumably this is what is intended; in which case it computes
"the unnormalized incomplete Gamma Function
 \Gamma(a,x) = \int_x^\infty dt t^{a-1} \exp(-t)
 for a real and x >= 0."

And again that is clear enough -- once you track it down!

In many places in the R documentation (including the "?" pages)
people have taken the trouble to spell out mathematical definitions
(where these can be given in reasonable space). Especially in cases
like the Incomplete Gamma and Beta functions, where there can be
dispute over what is meant (see above), it is surely wise to spell
it out!

Best wishes to all,
Ted.

>> On 31 Aug 2007, at 00:29, poolloopus at yahoo.com wrote:
>>
>>> Hello
>>>
>>> I am trying to evaluate an Incomplete gamma function
>>> in R. Library Zipfr gives the Igamma function. From
>>> Mathematica, I have:
>>>
>>> "Gamma[a, z] is the incomplete gamma function."
>>>
>>> In[16]: Gamma[9,11.1]
>>> Out[16]: 9000.5
>>>
>>> Trying the same in R, I get
>>>
>>>> Igamma(9,11.1)
>>> [1] 31319.5
>>> OR
>>>> Igamma(11.1,9)
>>> [1] 1300998
>>>
>>> I know I have to understand the theory and the math
>>> behind it rather than just ask for help, but while I
>>> am trying to do that (and only taking baby steps, I
>>> must admit), I was hoping someone could help me out.
>>>
>>> Regard
>>>
>>> Kris.

--------------------------------------------------------------------
E-Mail: (Ted Harding) <Ted.Harding at manchester.ac.uk>
Fax-to-email: +44 (0)870 094 0861
Date: 01-Sep-07                                       Time: 15:49:57
------------------------------ XFMail ------------------------------