[R] Regularized gamma function/ incomplete gamma function

[Martin Maechler]

    RV> I would be very grateful if you could help me with:

    RV> Given the regularized gamma function Reg=int_0^r
    RV> (x^(k-1)e^(-x))dx/int_0^Inf (x^(k-1)e^(-x))dx ; 0<r<Inf
    RV> (which is eventually the ratio of the Incomplete gamma
    RV> function by the gamma function),

and which is exactly what  R' s pgamma() is !

    RV> does anyone know of a
    RV> package in R that would evaluate the derivative of the
    RV> inverse of Reg with respect to k? I am aware that the
    RV> function "Rgamma.inv" of the package "Zipfr" evaluates
    RV> the inverse of Reg

[  well, the package's names is  'zipfR',   as you should know case
   does matter on decent computer environments

   Indeed, it seems that the author of zipfR has neither been aware
   that the (scaled / aka regularized) incomplete gamma (and beta,
   for that matter!) functions have been part of R all along.
   ...
   ... well , inspecting his code reveals he did know it.
   But why then on earth provide all the new <foo>gamma()
   functions, all trivially defined via pgamma(), qgamma() and
   gamma() ??
   Never mind ... Let's get to answer your Q ]
]

    RV> the inverse of Reg and I'm wondering wether there is a
    RV> function that would evaluate the derivative of the
    RV> inverse..

I'm a bit shocked by the lack of basic calculus knowledge both
in your question and even more in the answers.

I'm pretty sure that even before I've started studying math,
I knew the formula for getting the derivative of an inverse.

The mnemonic trick is  "dy / dx =   1/  (dx / dy)",
spelled out, that's

	d/dy f^{-1}(y) = 1 / f'(x) = 1 / f'(f^{-1}(y))

Now if you apply this to  f(x) = pgamma(x, a)
then the derivative of the inverse of the regularized incomplete
gamma function, i.e.e the derivative of qgamma()
is simply
	1 / dgamma(qgamma(x, a), a)

you can easily check this comparing with the results from
'numDeriv' if you want or just the simple one liner (computing
the difference ratio as approximate differential ratio) below:

For a = 1.25, and x = 0.2, e.g. :

  > sapply(10^-(3:9), function(e) diff(qgamma(.2 + c(-e,e), sh = 1.25))/(2*e))
  [1] 1.675105 1.675103 1.675103 1.675103 1.675103 1.675103 1.675103

  > 1/dgamma(qgamma(0.2, sh = 1.25), sh = 1.25)
  [1] 1.675103


Martin Maechler, ETH Zurich



    RV> Alternatively, a good numerical integration package/ or
    RV> simply a function that could evaluate the integral
    RV> int_0^r (log(x) x^(k-1) e^(-x))dx; 0<r<Inf would be
    RV> useful. I tried the function "int" of the package
    RV> "rmutil" but I'm not sure wether it is accurate
    RV> especially for small values of k. Does R have a powerful
    RV> numerical integration package that can deal with such
    RV> functions especially when the limit close to zero in +
    RV> or - Inf?

    RV> Many thanks for this opportunity to post our queries,

    RV> Amy