[R] Incomplete Gamma Functions and GammaDistribution Doc errata.

[Prof Brian Ripley]

Note corrections to two typos, since Master Custer seems to get confused
by tiny details.

On Fri, 19 Mar 2004, Prof Brian Ripley wrote:

> 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*.

the *shape*.

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

P(nu, x)  = pgamma(x, shape = nu) * P(nu, Inf)
          = pgamma(x, shape = nu) * 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