Inverse Gamma Distribution – How to Write the Derivative of the Inverse Gamma Function?

derivativeinverse gamma distributionmarkov-chain-montecarlo

I have recently been writing an R program on the inverse of the gamma function and the derivative of the inverse function. Now there is some confusion I would like to ask for advice.

I have written the program on how to find the inverse function of the gamma function.

invgamma=function(x,k)
{
  c=sqrt(2*pi)/exp(1)-gamma(k)
  L=log((x+c)/sqrt(2*pi))
  Wl=W(L/exp(1))
  igamma=L/Wl+1/2
  igamma
}

And how do I write the derivative of the inverse gamma? I used the following code, but I found that there are NaN in it

> invgade<-deriv(invgamma(x,2),"a")  
> invgade
expression({
    .value <- c(NaN, 4.28042788119302, 2.46198926331048, 3.58938901428173, 
    4.06875857308362, NaN, NaN, 2.99397817351358, 3.56516057177189, 
    2.18905646528529, 3.26296763434172, 3.26296763434172, 3.73900201030914, 
    2.56614173953932, 4.38384671521946, 2.46198926331048, 4.39838956859244, 
    2.18905646528529, 4.27174916443118, 3.87798639188068)
    .grad <- array(0, c(length(.value), 1L), list(NULL, c("a")))
    .grad[, "a"] <- 0
    attr(.value, "gradient") <- .grad
    .value
})

thank you very much!

Best Answer

If you have a function $f$ with inverse function $f^{-1}$, the derivative of $f^{-1}$ is:

$$ (f^{-1})'(x) = {1 \over f'(f^{-1}(x))}$$

from the chain rule. Note here $x$ is the argument to $f^{-1}$, not the argument to $f$. So, if you have the derivative of the gamma function handy, and the inverse gamma function, you can compute the derivative of the inverse gamma function.

The digamma function ($\Psi$) in R (statsmod package) computes the derivative of the log of the $\Gamma$ function, so evidently you can calculate $\Gamma '(z) = \Psi(z)\Gamma(z)$ with existing R functions. All you need to complete the job is something that calculates $\Gamma^{-1}$, which your function doesn't appear to do correctly (without W being defined.) You can either implement it using a truncated series representation, some other way, or just call the uniroot function and have it calculated numerically.

Sample code:

library(statsmod)

# For gamma(x) > 1
d_inverse_gamma <- function(x) {
  
  # function for root finding; refer to @whuber's
  # comment below for the rationale for using the log

  foo <- function(z) lgamma(z) - log(x)
  
  # Bracket root
  lb <- log(x)
  ub <- 2*lb
  while(foo(ub) < 0) {
    lb <- ub
    ub <- ub * 2
  }
  
  res <- uniroot(foo, lower=lb, upper=ub)
  digamma(res$root) *  x
}

which results in:

d_inverse_gamma(22.5)
[1] 33.67255

As a check, we can calculate the approximate derivative $(f(x+\delta)-f(x-\delta))/2\delta$ for a small delta. Repurposing the interior part of our function allows us to find the root, which is 4.957011, and we get:

(gamma(4.957011+0.001)-gamma(4.957011-0.001))/0.002
[1] 33.67247

very close to our answer above.

This latter result forces me to suggest that you just do it the easy (approximate derivative) way, if you don't need more than four or five decimal places of accuracy.

Related Question