R help - Mar 2004 - Incomplete Gamma Functions and GammaDistribution Doc errata.

Hello all,

In the course of trying to implement the CDF of an
InverseGammaDistribution, I have run across the need for an igamma()
function. Several others have needed this function but the answers I
have found so far are not totally clear to me. I'm writing for three
reasons:
	1) to present a small error in the docs
	2) to clarify the approach we are expected to take
	3) to request,for the ease of use of myself and others, that the 	  R
developers implement igamma(), document the approach they 	  want users
to take or do both.

1) An error in the docs of the gamma distribution:
-------------------------------------------------

The GammaDistribution(param1,param2) is defined alternatively with two
versions of the second parameter which are inverses of one another.
'Rate' arises from one derivation such as that presented in Mathworld: 
http://mathworld.wolfram.com/GammaDistribution.html
(see just below equation (2) )
        "define $\theta \equiv 1 / \lambda$ to be the time between
        changes"
which makes \theta in equation (3) a rate parameter. In the R
documentation (i.e. ?pgamma) this same equation is presented but the
parameter \theta is now called s and is termed the scale. This goes
against general usage of mathworld, wikipaedia, the COLT library, and
others who call the "s" term a "rate". 

Please either change the equation or replace the name with "rate" and
the label with r or nu or theta.



2) Using existing functions to build the igamma(\alpha,\beta):
-------------------------------------------------------------

The incomplete gamma has been mentioned in two previous emails that I
could find:
http://www.r-project.org/nocvs/mail/r-help/2000/1006.html
and
http://www.r-project.org/nocvs/mail/r-help/2000/3618.html
both of which suggest that the upper incomplete gamma can be obtained 
using the pgamma() function. 

Is Prof. Ripley's message (the first above) saying that

	upperIncompletegamma(x,\rate) = pgamma(x,\rate) * gamma(\rate) ?

I have not yet been able to see why this would be true. I keep tripping
myself up on the differing notation of R, mathworld, wikipaedia and
Abramowitz and Stegun's "Handbook of mathematical functions".

Would someone please confirm that this is mathematically correct? 


3) Adding an igamma() function or documentation for how to generate it.
--------------------------------------------------------------------

Since I am at least the third person not clever enough to see the way to
derive the igamma directly, because the lack of an igamma() function has
cost me a chunk of my day, because using functions related to
distributions to derive pure mathematical functions is counter
intuitive, and because two implementations already exist, I hope the R
developers would be willing to incorporate an igamma() function directly
into R. 

The implementation could either use the code posted to the list:
http://www.r-project.org/nocvs/mail/r-help/1998/0295.html
or be a convinience function around the pgamma and gamma solution
suggested in the emails presented above. 

Alternatively, if this approach doesn't suit you, would someone please
add an explanation to the documentation for the gamma function
explaining how to generate an igamma(\alpha, \beta) using the current
functions?

Thanks to all,
adrian

On Thu, 18 Mar 2004, Adrian Custer wrote:
> In the course of trying to implement the CDF of an
> InverseGammaDistribution, I have run across the need for an igamma()
> function. Several others have needed this function but the answers I
> have found so far are not totally clear to me. I'm writing for three
> reasons:
> 	1) to present a small error in the docs
> 	2) to clarify the approach we are expected to take
> 	3) to request,for the ease of use of myself and others, that the 	  R
> developers implement igamma(), document the approach they 	  want users
> to take or do both.
We would rather you defer to someone who leaps to fewer false assumptions.
> 1) An error in the docs of the gamma distribution:
> -------------------------------------------------
> 
> The GammaDistribution(param1,param2) is defined alternatively with two
> versions of the second parameter which are inverses of one another.
> 'Rate' arises from one derivation such as that presented in
Mathworld:
> http://mathworld.wolfram.com/GammaDistribution.html
> (see just below equation (2) )
>         "define $\theta \equiv 1 / \lambda$ to be the time between
>         changes"
Do please read what your reference says.  It calls \lambda the `rate of 
change', and that is the only occurrence of the word `rate' on the page.
> which makes \theta in equation (3) a rate parameter. In the R
That seems an impossible interpretation of what is actually said.
> documentation (i.e. ?pgamma) this same equation is presented but the
> parameter \theta is now called s and is termed the scale. This goes
> against general usage of mathworld, wikipaedia, the COLT library, and
> others who call the "s" term a "rate". 
`mathworld' plainly does not, and you have given us no others to verify.
> Please either change the equation or replace the name with "rate"
and
> the label with r or nu or theta.
Please look `rate' up in your dictionary, and look the gamma density up in 
a score of statistical references.  Your reference calls what R calls 
`rate' a `rate of change', and you insist it means that 1/rate is `a
rate
parameter'!

> 2) Using existing functions to build the igamma(\alpha,\beta):
> -------------------------------------------------------------
> 
> The incomplete gamma has been mentioned in two previous emails that I
> could find:
> http://www.r-project.org/nocvs/mail/r-help/2000/1006.html
> and
> http://www.r-project.org/nocvs/mail/r-help/2000/3618.html
> both of which suggest that the upper incomplete gamma can be obtained 
> using the pgamma() function. 
Neither mention the `upper incomplete gamma'.
> Is Prof. Ripley's message (the first above) saying that
> 
> 	upperIncompletegamma(x,\rate) = pgamma(x,\rate) * gamma(\rate) ?

Let me quote it to you, since by this point I have problems with your 
veracity.

  > Has anybody implemented the _incomplete_ gamma function in R?

  What do you think the cumulative distribution function of a gamma
  distribution is? (Now, there is a little room for disagreement here,
  but I think

  pgamma(x, nu)*gamma(nu)

  might be what you are looking for.) 

So

-  `Upper incomplete gamma' is not mentioned.
-  In your own reference, the parameters of an incomplete gamma do not 
contain a rate,

and so what makes you think nu is \rate?  The second argument of pgamma is
the *scale*.

I believe the `incomplete Gamma' of people like Pearson is what your 
reference calls the `lower incomplete Gamma', but note I did say

	there is a little room for disagreement here
and

	might be what you are looking for.

Let's look at Abramowitz and Stegun,

http://jove.prohosting.com/~skripty/page_260.htm

They have P(a, x) as the integral from 0 to x, whereas in Mathematica
notation \gamma(a, x) is the integral from 0 to a with shape parameter x.  
(That is pretty strange, since it is a that is normally varied for fixed
x.) So Abramowitz and Stegun's incomplete Gamma is the lower incomplete
Gamma of Mathematica with a more sensible set of parameter names.  Then in
your reference the gamma cdf for \lambda=1 is

D(x) = pgamma(x, shape = h) = P(h, x)/P(h, Inf)

and hence

P(nu, x)  = pgamma(x, shape = h) * P(h, Inf) 
          = pgamma(x, shape = h) * gamma(nu)

just as I suggested.
> I have not yet been able to see why this would be true. I keep tripping
> myself up on the differing notation of R, mathworld, wikipaedia and
> Abramowitz and Stegun's "Handbook of mathematical functions".
> 
> Would someone please confirm that this is mathematically correct? 
It is pretty discourteous to ask other people to confirm that an expert
posting is `mathematically correct', especially given the phrases you
selectively omitted.

> 3) Adding an igamma() function or documentation for how to generate it.
> --------------------------------------------------------------------
> 
> Since I am at least the third person not clever enough to see the way to
> derive the igamma directly, because the lack of an igamma() function has
> cost me a chunk of my day, because using functions related to
> distributions to derive pure mathematical functions is counter
> intuitive, and because two implementations already exist, I hope the R
> developers would be willing to incorporate an igamma() function directly
> into R. 
> 
> The implementation could either use the code posted to the list:
> http://www.r-project.org/nocvs/mail/r-help/1998/0295.html
> or be a convinience function around the pgamma and gamma solution
> suggested in the emails presented above. 
> 
> Alternatively, if this approach doesn't suit you, would someone please
> add an explanation to the documentation for the gamma function
> explaining how to generate an igamma(\alpha, \beta) using the current
> functions?
It is very easy, *once* you know what you mean by an incomplete Gamma
function.  Note that if you read mathworld, you are unlikely to be able to
communicate with either mathematicians or statisticians, even if you read
it accurately.

-- 
Brian D. Ripley,                  ripley at stats.ox.ac.uk
Professor of Applied Statistics,  http://www.stats.ox.ac.uk/~ripley/
University of Oxford,             Tel:  +44 1865 272861 (self)
1 South Parks Road,                     +44 1865 272866 (PA)
Oxford OX1 3TG, UK                Fax:  +44 1865 272595

>>>>> "Adrian" == Adrian Custer <acuster at
nature.berkeley.edu>
>>>>>     on Thu, 18 Mar 2004 18:47:54 -0800 writes:
    Adrian> Hello all,
    Adrian> In the course of trying to implement the CDF of an
    Adrian> InverseGammaDistribution, I have run across the need for an
igamma()
    Adrian> function. Several others have needed this function but the
answers I
    Adrian> have found so far are not totally clear to me. I'm writing
for three
    Adrian> reasons:

    Adrian> 1) to present a small error in the docs
I don't agree.
Did you look at the postscript or PDF version of the help file?
{you can get with
     help(pgamma, offline=TRUE)
 when you have a working latex installation}

If you compare its formula with Abramowitz & Stegun 6.5.1
its pretty pretty very very obvious

that   P(a,x) == pgamma(x,a)

{so forget about 'rate' or 'scale' which are both 1 by default!)

Why should this be hard to see? Shuldn't you first do your
homework before starting to insult the R developers?

Regards,
Martin Maechler, Seminar fuer Statistik, ETH Zurich Switzerland

My appologies to everyone. 

Prof. Ripley, please realize that I did not intend to suggest that your
mathematics were wrong. My question about mathematical correctness
refered to my presention of your formula as an equation. I have great
admiration for your work and its influences on my field and, therefore,
start with the presumption that all errors must be my own.

Mr. Maechler, please understand that I had no intent to "insult"
anyone.
I wanted to communicate three things:

        1) to point out what I understood to be an error in the
        documentation (which you both point out to be my own
        misunderstanding of the terminology).
        
        2) to ask for confirmation that I was interpreting past emails
        correctly.
        
        3) to argue for an R implementation of an incomplete gamma
        mathematical function.

To the R community, please accept my appologies for a message which you
may have taken to be insulting. I am deeply thankful for your work, for
your code and for this email list.

--adrian