Re: [R] Incomplete Gamma function
Hello Ted
thanks for the comments below.
You point out below some less-than-perfect aspects of
some of my documentation (bizarrely, this is not the first time that
this
has happened. The real problem is that *other people* insist on
reading the docs, when as everyone knows, the real purpose
of documentation is to remind the *developer* what he's done ;-).
The issue with the gsl_sf_ prefix is interesting. I decided that
having an R function named "gsl_sf_gamma_inc()" would be
a bit long-winded, so I removed the prefix. I did this elsewhere
too, notably Bessel.Rd, because I felt that the GSL naming scheme
was sufficiently distinct from R's in the case of Bessel functions to
obviate typing "gsl_sf_" every single time.
But as you point out, this is potentially confusing to a user; as a
developer, one tends to see the Rd file as one of a number of files
that one must edit, and in the case of the gsl library, the R and C
code is very regimented and has adequate inline comments. But
I'd lost sight of the fact that many users' sole source of documentation
is the Rd files. And in the case of gsl, the documentation is just a
pointer to ams-55 via gsl-ref [the GSL reference manual]. I now
suspect that my short function definitions have unnecessarily
obscured that audit trail.
What I think I'll do (read "will add to my to-do-list")
is to (eg) define R function gsl_sf_gamma_inc()
and then define gamma_inc <- gsl_sf_gamma_inc
and add appropriate aliases to the docs.
Then I could deprecate gsl_inc(); and then possibly defunct them.
Anyone got any comments on this? Do other package developers
have any experiences of making functions defunct that would be
interesting?
How long do folk leave functions deprecated before defuncting them?
best wishes
Robin
On 1 Sep 2007, at 15:50, (Ted Harding) wrote:
[snip]
>
> 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, [EMAIL PROTECTED] 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) <[EMAIL PROTECTED]>
> Fax-to-email: +44 (0)870 094 0861
> Date: 01-Sep-07 Time: 15:49:57
> -- XFMail --
>
> __
> R-help@stat.math.ethz.ch mailing list
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide http://www.R-project.org/posting-
> guide.html
> and provide commented, minimal, self-contained, reproducible code.
--
Robin Hankin
Uncertainty Analyst
National Oceanography Centre, Southampton
European Way, Southampton SO14 3ZH, UK
tel 023-8059-7743
__
R-help@stat.math.ethz.ch mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
Re: [R] Incomplete Gamma function
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, [EMAIL PROTECTED] wrote:
>>
>>> Hello
>>>
>>> I am trying to evaluate an Incomplete gamma function
>>> in R. Library Zipfr gives the Igamma function. From
>>> Mathematica, I have:
>>>
>>> "Ga
Re: [R] Incomplete Gamma function
> "AN" == Anup Nandialath <[EMAIL PROTECTED]> > on Fri, 31 Aug 2007 13:15:08 -0700 (PDT) writes: AN> Hi Kris, You just need to understand the mathematics of AN> the incomplete gamma function and the various AN> relationships it has. The answers from both Mathematica AN> and R are correct, except that they are giving you AN> different estimated quantities. It depends on the way AN> the gamma function is written. For instance in R to get AN> the equivalent result from mathematica you should do the AN> following AN> answer <- gamma(9) - Igamma(9,11.1). This will give you AN> the incomplete gamma for (9,11.1) as given by AN> Mathematica. AN> You can read more about the model and am sure you will AN> figure it out. Yes, and then, *PLEASE* do as Brian Ripley suggested, and understand that the (normalized) incomplete gamma function is the same as the gamma distribution functions and has been available in S and R "for ever" __and__ in R has been very thoroughly tested and quite a few extreme cases have been made to work more accurately, etc etc: pgamma() ! BTW: The same applies to the incomplete beta function which -- in one of it's equivalent forms -- is called the beta distribution function in probability and statistics and has been available in S and R "for ever", and for R, has been very carefully tested and for extreme border cases gradually improved over the years, most recently for the upcoming R 2.6.0, where the precision of pbeta(*, log=TRUE) has been dramatically improved in some extreme tail/range cases. --> which benefits pt(), pf(), pbinom() etc in equivalent situations. Martin Maechler, ETH Zurich and R-core AN> Anup AN> - Got a little couch AN> potato? Check out fun summer activities for kids. AN> [[alternative HTML version deleted]] AN> __ AN> R-help@stat.math.ethz.ch mailing list AN> https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do AN> read the posting guide AN> http://www.R-project.org/posting-guide.html and provide AN> commented, minimal, self-contained, reproducible code. __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
Re: [R] Incomplete Gamma function
Hi Kris, You just need to understand the mathematics of the incomplete gamma function and the various relationships it has. The answers from both Mathematica and R are correct, except that they are giving you different estimated quantities. It depends on the way the gamma function is written. For instance in R to get the equivalent result from mathematica you should do the following answer <- gamma(9) - Igamma(9,11.1). This will give you the incomplete gamma for (9,11.1) as given by Mathematica. You can read more about the model and am sure you will figure it out. Regards Anup - Got a little couch potato? Check out fun summer activities for kids. [[alternative HTML version deleted]] __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
Re: [R] Incomplete Gamma function
On 31 Aug 2007, at 14:06, 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. GSL gives both tails for the normalized incomplete gamma function but not the unnormalized version, for which only the upper incomplete gamma function is given. This seems a strange omission; I'll raise it on the GSL mailing list. re your request for clarification in the gsl helpfiles: I'm uncertain about this. The gsl library is supposed to be a "no-brainer" wrapper for GSL functions; anything documented in .Rd files will be 100% dominated by gsl-ref [the GSL manual]. That's why the documentation in gsl is restricted to the rather curt "see the GSL documentation". But I appreciate that this seems rather unhelpful. Best wishes rksh -- Robin Hankin Uncertainty Analyst National Oceanography Centre, Southampton European Way, Southampton SO14 3ZH, UK tel 023-8059-7743 __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
Re: [R] Incomplete Gamma function
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. > On 31 Aug 2007, at 00:29, [EMAIL PROTECTED] 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. >> >> >> >> __ >> __ >> Got a little couch potato? >> Check out fun summer activities for kids. >> >> __ >> R-help@stat.math.ethz.ch mailing list >> https://stat.ethz.ch/mailman/listinfo/r-help >> PLEASE do read the posting guide http://www.R-project.org/posting- >> guide.html >> and provide commented, minimal, self-contained, reproducible code. > > -- > Robin Hankin > Uncertainty Analyst > National Oceanography Centre, Southampton > European Way, Southampton SO14 3ZH, UK > tel 023-8059-7743 > > __ > R-help@stat.math.ethz.ch mailing list > https://stat.ethz.ch/mailman/listinfo/r-help > PLEASE do read the posting guide http://www.R-project.org/posting-guide.html > and provide commented, minimal, self-contained, reproducible code. > -- Brian D. Ripley, [EMAIL PROTECTED] 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, UKFax: +44 1865 272595 __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
Re: [R] Incomplete Gamma function
Hi Kris lgamma() gives the log of the gamma function. You need gamma_inc() of the gsl package, a wrapper for the GSL library: > gamma_inc(9,11.1) [1] 9000.501 > HTH rksh On 31 Aug 2007, at 00:29, [EMAIL PROTECTED] 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. > > > > __ > __ > Got a little couch potato? > Check out fun summer activities for kids. > > __ > R-help@stat.math.ethz.ch mailing list > https://stat.ethz.ch/mailman/listinfo/r-help > PLEASE do read the posting guide http://www.R-project.org/posting- > guide.html > and provide commented, minimal, self-contained, reproducible code. -- Robin Hankin Uncertainty Analyst National Oceanography Centre, Southampton European Way, Southampton SO14 3ZH, UK tel 023-8059-7743 __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
[R] Incomplete Gamma function
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. Got a little couch potato? Check out fun summer activities for kids. __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
Re: [R] Incomplete gamma function (was "integrate()", {was "mathematica -> r ..."})
> "MM" == Martin Maechler <[EMAIL PROTECTED]>
> on Tue, 8 Aug 2006 09:55:50 +0200 writes:
> "Leo" == Leo Gürtler <[EMAIL PROTECTED]>
> on Tue, 08 Aug 2006 00:13:19 +0200 writes:
Leo> Dear R-list,
Leo> I try to transform a mathematica script to R.
Leo> ###relevant part of the Mathematica script
Leo> (* p_sv *)
Leo> dd = NN (DsD - DD^2);
Leo> lownum = NN (L-DD)^2;
Leo> upnum = NN (H-DD)^2;
Leo> low = lownum/(2s^2);
Leo> up = upnum/(2s^2);
Leo> psv = NIntegrate[1/(s^NN) Exp[-dd/(2s^2)]
Leo>(Gamma[1/2,0,up] + Gamma[1/2,0,low]),{s,sL,sH}, MinRecursion-> 3];
Leo>PSV = psv/Sqrt[2NN];
Leo> Print["- Results "];
Leo> Print[" "];
Leo> Print["p(sv|D_1D_2I) = const. ",N[PSV,6]];
Leo>
Leo> # R part
Leo> library(fOptions)
Leo> ###raw values for reproduction
Leo> NN <- 58
Leo> dd <- 0.411769
Leo> lownum <- 20.81512
Leo> upnum <- 6.741643
Leo> sL <- 0.029
Leo> sH <- 0.092
Leo> ###
Leo> integpsv <- function(s) { 1 / (s^NN) * exp(-dd / (2 * s^2)) *
Leo>( (igamma((upnum/(2*s^2)),1/2) - igamma(0,1/2) ) +
Leo>(igamma((lownum/(2*s^2)),1/2) - igamma(0,1/2) ) )
Leo> }
Leo> psv <- integrate(integpsv, lower=sL, upper=sH)
Leo> PSV <- psv$value / sqrt(2*NN)
Leo> print("- Results \n")
Leo> print(paste("p(sv|D_1D_2I) = const. ",PSV, sep=""))
Leo> The results of variable "PSV" are not the same.
Leo> In mathematica -> PSV ~ 2.67223e+47
Leo> with rounding errors due to the initial values, in R -> PSV ~ 1.5e+47
Leo> I am not that familiar with gamma functions and integration, thus I
Leo> assume there the source of the problem can be located.
MM> Yes.
MM> A few remarks
MM> 1) No need to use package "fOptions" and igamma();
MM> standard R's pgamma() is all you need
MM> {igamma() has added value only for *complex* arguments!}
MM> 2) igamma(0, 1/2) == pgamma(0, 1/2) == 0 , so you can really
MM> drop them from your integrand.
MM> integpsv <- function(s) {
MM>1 / (s^NN) * exp(-dd / (2 * s^2)) *
MM>( pgamma(upnum/(2*s^2), 1/2) + pgamma(lownum/(2*s^2), 1/2) )
MM> }
[]
MM> However, if I experiment, using integrate() in two parts, or using many
other
MM> numerical integration approximators,
MM> all methods give ( your 'psv', not PSV )
MM> integrate(integpsv, lower=sL, upper=sH)
MM> a value of 1.623779e+48 (which leads to your PSV of 1.5076e+47)
MM> Could it be that you are not using the same definition of
MM> incomplete gamma in Mathematica and R ?
Offlist, Leo sent me Mathematica's definition
of
Gamma[a, z0, z1] := integral_z0^z1 t^(a-1) exp(-t) dt
Now if you compare this with what ?pgamma (not ?gamma !) tells you,
namely that R uses (Abramowitz and Stegun's definition of the
incomplete gamma function)
pgamma(x, a) = 1/ Gamma(a) * integral_0^x t^(a-1) exp(-t) dt
which has a normalizing factor: In your case above, it is
Gamma(1/2) with its well-known value of sqrt(pi).
And indeed, if you multiply the current result by sqrt(pi), you
get what you want -- and did get from Mathematica:
> (1.623779e+48 / sqrt(2*NN)) * sqrt(pi)
[1] 2.672224e+47
Regards,
Martin Maechler, ETH Zurich
__
R-help@stat.math.ethz.ch mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
