如果不进行并行化，我们可以通过找到所需函数的高性能版本来加快速度。dqrng package的runif版本比base快3倍，在长向量上sum(x) / length(x)比mean(x)快一点。

library(dqrng)
nn = 1e6
bench::mark(
  mean(runif(nn, 1, 5)),
  mean(dqrunif(nn, 1, 5)),
  sum(dqrunif(nn, 1, 5)) / nn,
  check = FALSE
)
# # A tibble: 3 × 13
#   expression                     min   median `itr/sec` mem_alloc `gc/sec` n_itr  n_gc
#   <bch:expr>                <bch:tm> <bch:tm>     <dbl> <bch:byt>    <dbl> <int> <dbl>
# 1 mean(runif(nn, 1, 5))      37.09ms     40ms      25.1    7.63MB     5.02    10     2
# 2 mean(dqrunif(nn, 1, 5))     8.23ms     12ms      85.2    7.63MB    11.4     30     4
# 3 sum(dqrunif(nn, 1, 5))/nn   7.78ms    8.5ms     117.     7.63MB    13.3     44     5
# # ℹ 5 more variables: total_time <bch:tm>, result <list>, memory <list>, time <list>,
# #   gc <list>

结合这两种方法，我们得到了> 3倍的速度提升，如果结合并行化，速度会更快，如neilfws的回答：

library(doParallel)
registerDoParallel(7)
system.time( { times(n) %dopar% sum(dqrunif(nn, 1, 5)) / nn } )
#    user  system elapsed 
# 178.940  31.608  86.432

不到1.5分钟（相比之下，我的笔记本电脑上的原始代码大约需要6.5分钟）。

在R中有很多方法可以进行并行计算。您可以查看：

• Parallel Computation
• CRAN Task View: High-Performance and Parallel Computing with R

举个例子，在我的工作机器（一台普通的Surface Book 2）上使用doParallel库：

library(doParallel)

registerDoParallel(7)

# original version
system.time ( { replicate(10000, mean(runif(1e+06,1,5))) } )
 user  system elapsed 
 319.70   20.36  340.39

# parallel version 7 cores
 system.time( { times(10000) %dopar% mean(runif(1e+06,1,5)) } )
  user  system elapsed 
  6.06    1.14  125.75

所以大约2分钟，而不是5分钟多一点（不完全是“永远”，但足够长）。
其中一些也可能有所帮助。

贝茨分布可以看作是正态分布的多项式近似（当n接近Inf时，它是正态分布）。对于n = 1e6，正态分布是一个非常好的近似。使用1e5示例演示：

library(parallel)

# Direct computation of Bates r.v.
cl <- makeCluster(parallel::detectCores() - 1L, type = "PSOCK")
clusterEvalQ(cl, library(dqrng))
system.time(x1 <- unlist(parLapply(cl, 1:1e3, \(i) replicate(100, sum(dqrunif(1e6, 1, 5)))))/1e6)
#>    user  system elapsed 
#>    0.02    0.00  103.25

# normal approximation
x2 <- rnorm(1e5, 3, sqrt(16/12/1e6))

# Kolmogorov-Smirnov test for a difference between the two distributions
ks.test(x1, x2)
#> 
#>  Asymptotic two-sample Kolmogorov-Smirnov test
#> 
#> data:  x1 and x2
#> D = 0.00312, p-value = 0.7151
#> alternative hypothesis: two-sided

# plot the empirical CDFs
plot(ecdf(x1), col = "blue")
plot(ecdf(x2), col = "orange", add = TRUE)

这两组样本在数值上基本上是不可区分的。为了进行比较，绘制正态分布中两个不同样本的经验CDF。

plot(ecdf(rnorm(1e5, 3, sqrt(4/3/1e6))), col = "blue")
plot(ecdf(rnorm(1e5, 3, sqrt(4/3/1e6))), col = "orange", add = TRUE)

你不能在这里模拟贝茨分布。整数参数为n的Bates分布是（0，1）上n个均匀随机变量的均值。一个快速获得它的方法是：

n <- 6
nsims <- 100000
usims <- matrix(runif(n*nsims), nrow = n, ncol = nsims)
bates_sims <- colMeans(usims)