Inverse of gamma function with R

Question

With gamma() in R you get output 24 from input 5. How can I make inverse function which can get output 5 from input 24?

>gamma(5)
24

I found [invgamma] package, but I can’t understand if it is appropriate one nor how to use it…

Diffract

The invgamma package is for the inverse gamma distribution, rather than inverting the gamma function. – Undermine

Answer 1

Use uniroot to find the root of gamma(x) - a, where a is your number.

invgamma <- function(a, tol = .Machine$double.eps^0.5) {
  uniroot(\(x) gamma(x) - a, c(1, 1e2), extendInt = "upX", tol = tol)
}

invgamma(24)$root
#> [1] 5

invgamma(24)
#> $root
#> [1] 5
#> 
#> $f.root
#> [1] -5.207212e-11
#> 
#> $iter
#> [1] 19
#> 
#> $init.it
#> [1] NA
#> 
#> $estim.prec
#> [1] 7.450582e-09

^{Created on 2024-09-29 with reprex v2.1.0}

You can vectorize this function.

igamma <- Vectorize(invgamma, "a")
igamma(c(6, 24, 120)) |> t()
#>      root f.root        iter init.it estim.prec  
#> [1,] 4    -2.771117e-13 17   NA      7.450582e-09
#> [2,] 5    -5.207212e-11 19   NA      7.450582e-09
#> [3,] 6    2.727063e-11  18   NA      7.450583e-09

^{Created on 2024-09-29 with reprex v2.1.0}

Neoclassic

Answer 2

Below are two fully vectorized options using Newton-Raphson. The first is in base R. The second is a bit faster but requires the lamW package. Both take less than a second to calculate 1M inverses.

Depending on the usage, you may want additional checks, such as handling NA or nulls.

Note that x = gamma(y) has a minimum (x = 0.885603194410889) at y = 1.46163214496836 and approaches Inf as y approaches 0 and Inf, so there are two positive solutions to x = gamma(y) for x > 0.885603194410889. The functions below provide the value of y for y > 1.46163214496836 (the "principal branch").

---

Base R

This Q&A provides a decent initial approximation to begin the Newton-Raphson iterations.

gammainv <- function(x) {
  i <- 1:length(x)
  out <- numeric(length(x))
  out[b <- (x < 0.885603194410889)] <- NA # minimum of gamma function
  i <- i[!b]
  x <- x[i]
  out[i[b <- (x == 1)]] <- 2 # principal branch
  i <- i[!b]
  x <- x[!b]
  x <- log(x)
  # initial estimate
  m <- x/log1p(x) + 1
  k <- (lgamma(m) - x)/log(m)
  y <- m - k - 1 + (x - lgamma(m - k))/log(m - k)
  y[y < 1.46163214496836] <- 1.5 # principal branch
  # Newton-Raphson
  while (length(i)) {
    y0 <- y
    y <- y - (lgamma(y) - x)/digamma(y)
    b <- which(abs(y - y0) < 1e-13)
    if (length(b)) {
      out[i[b]] <- y[b]
      y <- y[-b]
      i <- i[-b]
      x <- x[-b]
    }
  }
  out
}

Testing:

1e6 - gamma(x <- gammainv(1e6))
#> [1] 2.211891e-09
dput(gamma(x))
#> 999999.999999998

gammainv(cumprod(1:170))
#>   [1]   2   3   4   5   6   7   8   9  10  11  12  13  14  15  16  17  18  19
#>  [19]  20  21  22  23  24  25  26  27  28  29  30  31  32  33  34  35  36  37
#>  [37]  38  39  40  41  42  43  44  45  46  47  48  49  50  51  52  53  54  55
#>  [55]  56  57  58  59  60  61  62  63  64  65  66  67  68  69  70  71  72  73
#>  [73]  74  75  76  77  78  79  80  81  82  83  84  85  86  87  88  89  90  91
#>  [91]  92  93  94  95  96  97  98  99 100 101 102 103 104 105 106 107 108 109
#> [109] 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127
#> [127] 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145
#> [145] 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160 161 162 163
#> [163] 164 165 166 167 168 169 170 171

x <- exp(runif(1e6, log(0.885603194410889), 709))
system.time(y <- gammainv(x))
#>    user  system elapsed 
#>    0.85    0.03    0.87
all.equal(x, gamma(y))
#> [1] TRUE

---

lamW

Additional speed using the Lambert W function, which provides a better initial approximation.

library(lamW)

gammainv <- function(x) {
  i <- 1:length(x)
  out <- numeric(length(x))
  out[b <- (x < 0.885603194410889)] <- NA # minimum of gamma function
  i <- i[!b]
  logsqrt2pi <- log(2*pi)/2
  x <- pmax(x, exp(logsqrt2pi - 1))
  x <- log(x[i])
  y <- exp(lambertW0((x - logsqrt2pi)/exp(1)) + 1) + 0.5
  # Newton-Raphson
  while (length(i)) {
    y0 <- y
    y <- y - (lgamma(y) - x)/digamma(y)
    b <- which(abs(y - y0) < 1e-13)
    if (length(b)) {
      out[i[b]] <- y[b]
      y <- y[-b]
      i <- i[-b]
      x <- x[-b]
    }
  }
  out
}

Testing

1e6 - gamma(y <- gammainv(1e6))
#> [1] 4.074536e-09
dput(gamma(y))
#> 999999.999999996

gammainv(cumprod(1:170))
#>   [1]   2   3   4   5   6   7   8   9  10  11  12  13  14  15  16  17  18  19
#>  [19]  20  21  22  23  24  25  26  27  28  29  30  31  32  33  34  35  36  37
#>  [37]  38  39  40  41  42  43  44  45  46  47  48  49  50  51  52  53  54  55
#>  [55]  56  57  58  59  60  61  62  63  64  65  66  67  68  69  70  71  72  73
#>  [73]  74  75  76  77  78  79  80  81  82  83  84  85  86  87  88  89  90  91
#>  [91]  92  93  94  95  96  97  98  99 100 101 102 103 104 105 106 107 108 109
#> [109] 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127
#> [127] 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145
#> [145] 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160 161 162 163
#> [163] 164 165 166 167 168 169 170 171

x <- exp(runif(1e6, log(0.885603194410889), 709))
system.time(y <- gammainv(x))
#>    user  system elapsed 
#>    0.73    0.05    0.56
all.equal(x, gamma(y))
#> [1] TRUE

Undermine

interesting and efficient implementation of inverse Gamma, +1!. Do you have any algorithmic reference about the code? Better to show it in the solution in case other wonder about how/why it works. – Inconstant

I added a couple links that discuss the initial approximations. – Undermine

That looks great, thanks for adding the references. – Inconstant