创建返回伽马分布的函数的更好方法？

Question

所以我创建了一个函数来生成伽马分布的数据帧。

目前，我有。

sample_gamma <- function(alpha,beta,n,iter) {
gamma.df <- as.data.frame(matrix(nrow = iter, ncol = 3))
colnames(gamma.df) <- c("iteration","mean","standard dev")
gamma.df$iteration <- c(1:iter)
  for (i in 1:iter) {
    gamma.dist <- rgamma(n,shape = alpha, rate = beta, scale = 1/beta)
    gamma.df[i,2] <- mean(gamma.dist)
    gamma.df[i,3] <- sd(gamma.dist)
  }
print(gamma.df)
}

该功能可以完成我需要做的所有事情，但我想知道是否有任何替代或更清洁的方法可以做到这一点

Answer 1

我将创建一个函数，该函数在一次迭代中返回 mean 和 sd。

sample_gamma <- function(alpha,beta,n) {
    dist <- rgamma(n,shape = alpha, rate = beta)
    c(mean = mean(dist), sd = sd(dist))
}

然后使用replicate重复它

t(replicate(5, sample_gamma(2, 3, 4)))

#          mean        sd
#[1,] 0.5990206 0.2404226
#[2,] 0.6108976 0.3083426
#[3,] 1.0616542 0.4602403
#[4,] 0.3415355 0.1543885
#[5,] 1.0558066 0.9659599

Answer 2

虽然我认为 Ronak Shah 的答案很简单且相对惯用（R-wise），但当扩展到高 iter 计数时，这里有一个更有效的答案（因为它只执行一次随机拉取）：

sample_gamma <- function(alpha, beta, n, iter) {
  mtx <- matrix(rgamma(n*iter, shape=alpha, rate=beta), nrow=n, ncol=iter)
  t(apply(mtx, 2, function(a) c(mean=mean(a), sd=sd(a))))
}
sample_gamma(2, 3, 4, 5)
#           mean         sd
# [1,] 0.6486220 0.22900833
# [2,] 0.8551055 0.07874287
# [3,] 0.7854750 0.72694260
# [4,] 0.7045878 0.24834502
# [5,] 1.1783301 0.25210538

基准测试：

microbenchmark::microbenchmark(
  RS=t(replicate(5, sample_gamma_RS(2,3,4))),
  r2=sample_gamma_r2(2,3,4,5)
)
# Unit: microseconds
#  expr   min     lq    mean median    uq    max neval
#    RS 413.7 493.70 757.884 743.80 946.1 1611.6   100
#    r2 405.2 461.15 681.630 706.35 898.6 1348.2   100

microbenchmark::microbenchmark(
  RS=t(replicate(500, sample_gamma_RS(2,3,4))),
  r2=sample_gamma_r2(2,3,4,500)
)
# Unit: milliseconds
#  expr    min       lq     mean   median       uq      max neval
#    RS 31.271 40.58735 56.44298 57.85735 65.08605  95.1866   100
#    r2 29.110 38.81230 53.99426 57.45820 61.35720 100.5820   100

microbenchmark::microbenchmark(
  RS=t(replicate(500, sample_gamma_RS(2,3,400))),
  r2=sample_gamma_r2(2,3,400,500)
)
# Unit: milliseconds
#  expr     min       lq     mean   median       uq      max neval
#    RS 60.6782 101.3112 121.3533 116.7464 140.8845 227.1904   100
#    r2 66.3892  81.0329 106.9920  98.7170 126.7742 198.3947   100

我承认我认为这会在性能上产生更大的差异。