答案here是优雅和快速的，但如果你不关心效率，

my_rpois <- function(n, lambda) {
   qpois(runif(n), lambda)
}

应该可以实现逆CDF方法。（lambda可以是向量）

修改我在你链接的帖子中的回答，以获取一个率向量：

r.poisson1 <- function(n, l = 0.5) {
  U <- runif(n)
  i <- 0L
  X <- integer(n)
  p <- exp(-l)
  F <- p
  idx <- 1:n
  while (length(idx)) {
    bln <- U < F
    X[idx[bln]] <- i
    idx <- idx[!bln]
    U <- U[!bln]
    l <- l[!bln]
    p <- p[!bln]
    F <- F[!bln]
    p <- l*p/(i <- i + 1L)
    F <- F + p
  }
  X
}

测试：

n <- 10
beta0 <- 1
beta1 <- 0.2
x <- runif(n=n, min=0, max=1.5)
mu <- exp(beta0 + beta1 * x)
system.time({y <- replicate(1e5, r.poisson1(n, mu))})
#>    user  system elapsed 
#>    3.30    0.03    3.49
system.time({z <- replicate(1e5, rpois(n, mu))})
#>    user  system elapsed 
#>    0.45    0.03    0.50
data.frame(
  mean.y = rowMeans(y),
  mean.z = rowMeans(z),
  var.y = apply(y, 1, var),
  var.z = apply(z, 1, var)
)
#>     mean.y  mean.z    var.y    var.z
#> 1  3.06219 3.05722 3.064813 3.067737
#> 2  2.73974 2.73555 2.748812 2.734944
#> 3  2.93172 2.93107 2.965167 2.913688
#> 4  3.03028 3.03886 3.031133 3.043200
#> 5  3.43201 3.43630 3.444752 3.447537
#> 6  3.47632 3.46268 3.488834 3.436982
#> 7  3.08662 3.09912 3.069428 3.108086
#> 8  3.16495 3.16348 3.169913 3.152986
#> 9  2.74217 2.74010 2.738821 2.754580
#> 10 2.82742 2.83341 2.813164 2.836006