以下是两个使用 Newton-Raphson 的完全矢量化选项。第一个是在基础 R 中。第二个稍快一些，但需要包lamW。两者都花费不到一秒钟来计算 1M 逆。

根据使用情况，您可能需要进行额外的检查，例如处理NA或空值。

请注意，在 处x = gamma(y)有最小值（x = 0.885603194410889） ，并且当接近和时y = 1.46163214496836接近，因此 有两个正解。下面的函数提供了 的值（“主分支”）。Infy0Infx = gamma(y)x > 0.885603194410889yy > 1.46163214496836

基R

为开始牛顿-拉夫森迭代提供了一个不错的初始近似值。

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
}

测试：

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

使用 Lambert W 函数来提高速度，从而提供更好的。

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
}

测试

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