如何在没有循环的情况下进行模拟？

Posted 2023-03-12

【中文标题】如何在没有循环的情况下进行模拟？

【英文标题】：How to do simulation without loop?

【发布时间】：2021-03-31 12:25:24

【问题描述】：

我正在编写一个模拟函数来计算R 中的 t 检验的功效。但是，在R 中编写循环效率不高，有没有其他方法可以在没有循环的情况下实现我的目标？

#Define a simulation function
simulation <- function(N,alpha,sigma,diff,mu1)
  p_values = c()
  for(i in 1:10000)
    group1 <- rnorm(n=N/2, mean = mu1, sd = sigma)
    group2 <- rnorm(n=N/2, mean = mu1+diff, sd = sigma)
    p_values[i] <- t.test(group1,group2)$p.value
  
  
  prop.table(table(p_values<alpha))[["TRUE"]]

【问题讨论】：

我相信循环在 R 中不一定很慢。使用“lapply”变体也不一定更快。与矢量化相比，循环很慢！如果没有很好的方法来矢量化您的问题，那么选择循环或“lapply”解决方案（这也只是对作为参数传递的矢量/列表中元素的循环非常好。

***.com/questions/28983292/…

感谢@BenBolker 的那个兔子洞。所以最后似乎都归结为哪些部分是用 C 有效编写的，哪些部分不是的细节？我倾向于自己编写 lapply 解决方案，因为它们不会使代码混乱，并且功能可以更容易隐藏/重用。但是反身的“不要在 R 中编写循环”似乎有些言过其实。

【参考方案1】：

tl;dr 循环很好。我发现显着加快这一速度的唯一方法是编写一个定制版本，将stats:::t.test.default 剥离为仅计算 p 值所需的基本代码（跳过不同选项的测试、计算置信区间等）。这获得了大约 2 倍的加速；如果不使用 C++ 编码（例如使用 Rcpp 包），我没有看到一种简单的方法来进一步加快速度。

更多注释：

预分配p_values 向量是我尝试的第一件事，但总体差别不大（t.test() 函数是瓶颈） 将prop_table(...) 替换为mean(p_values<alpha) 也没什么不同 power.t.test() 解决了同样的问题（我认为：

我不确定它是否假设方差相等

）并且速度快得多（但这可能不是重点你的问题） 另一种加快速度的可能方法（尽管我怀疑它会做很多事情）是一次选择所有正常偏差并将它们粘贴到适当尺寸的矩阵中，然后索引矩阵（而不是每次都调用rnorm() ）。这看起来很烦人，我猜在这种情况下不会有太大的不同。 您实际上可以对整个计算进行矢量化 - 但如果您想做一些更专业/更复杂的事情，这可能行不通。我写了它：它没有给出与其他模拟人生完全相同的答案（0.3498 vs 0.0.35215 for nsim=1e5），但我认为呢？这是因为随机数的分配顺序略有不同。令我惊讶的是，它并没有比第二版快多少……

## rewrite sim1 function slightly: add convenient default values
sim1 <- function(N=1000,alpha=0.05,sigma=1,diff=0.1,mu1=1, nsim=1e4) 
    p_values = c()
    set.seed(12345)
    for(i in seq(nsim)) 
        group1 <- rnorm(n=N/2, mean = mu1, sd = sigma)
        group2 <- rnorm(n=N/2, mean = mu1+diff, sd = sigma)
        p_values[i] <- t.test(group1,group2)$p.value
    
    prop.table(table(p_values<alpha))[["TRUE"]]

## stripped-down function to compute t-test p value, based on stats:::t.test.default
my_t <- function(x,y) 
    vx <- var(x)
    vy <- var(y)
    nx <- length(x)
    ny <- length(y)
    stderrx <- sqrt(vx/nx)
    stderry <- sqrt(vy/ny)
    stderr <- sqrt(stderrx^2 + stderry^2)
    df <- stderr^4/(stderrx^4/(nx - 1) + stderry^4/(ny - 1))
    tstat <- (mean(x) - mean(y))/stderr
    pval <- 2 * pt(-abs(tstat), df)
    return(pval)

## faster sim function
sim2 <- function(N=1000,alpha=0.05,sigma=1,diff=0.1,mu1=1, nsim=1e4) 
    p_values <- numeric(nsim)  ## pre-allocate loop
    set.seed(12345)
    for(i in seq(nsim)) 
        group1 <- rnorm(n=N/2, mean = mu1, sd = sigma)
        group2 <- rnorm(n=N/2, mean = mu1+diff, sd = sigma)
        p_values[i] <- my_t(group1, group2)
    
    mean(p_values<alpha) ## replace prop.table with cheaper alternative

## vectorized sim function
sim3 <- function(N=1000,alpha=0.05,sigma=1,diff=0.1,mu1=1, nsim=1e4) 
    set.seed(12345)
    group1 <- matrix(rnorm(n=N/2*nsim, mean=mu1,sd=sigma),
                     nrow=nsim)
    group2 <- matrix(rnorm(n=N/2*nsim, mean=mu1+diff,sd=sigma),
                     nrow=nsim)
    vx <- apply(group1, 1, var)
    vy <- apply(group2, 1, var)
    nx <- ny <- N/2
    stderrx <- sqrt(vx/nx)
    stderry <- sqrt(vy/ny)
    stderr <- sqrt(stderrx^2 + stderry^2)
    df <- stderr^4/(stderrx^4/(nx - 1) + stderry^4/(ny - 1))
    tstat <- (rowMeans(group1) - rowMeans(group2))/stderr
    p_values <- 2 * pt(-abs(tstat), df)
    mean(p_values<alpha) ## replace prop.table with cheaper alternative


identical(sim1(), sim2())  ## TRUE (value= 0.3553)
system.time(sim1(nsim=1e5))  ## 11.6 seconds
system.time(sim2(nsim=1e5))  ## 6 seconds
power.t.test(n=500,delta=0.1,sd=1)  ## value=0.3518

【参考方案2】：

您实际上可以完全退出 for 循环，方法是生成 group1 和 group2 的单个流，然后将它们标注为 ncol = 10000 的矩阵，然后 sapply 处理 t.test：

bigNum <- 100000
iters <- bigNum * (N/2)
group1 <- rnorm(n=iters, mean = mu1, sd = sigma)
group2 <- rnorm(n=iters, mean = mu1+diff, sd = sigma)
m1 <- matrix(group1, ncol = bigNum)
m2 <- matrix(group2, ncol = bigNum)
pvalues <- sapply(1:bigNum, function(x) t.test(m1[ , x], m2[ , x])$p.value)

【讨论】：

是的，但我认为这实际上不会加快速度。您尝试过计时测试吗...？

没有测试。但是您可以先在 bigNum = 1000 处尝试。让我们知道哪种解决方案适合您（如果有）。祝你好运。

我得到这个解决方案比我系统上的原始示例花费 longer（bigNum=1e5 为 13+ vs 11+ 秒），可能是因为内存分配花费了时间... 而我的两个解决方案大约需要 6 秒

【参考方案3】：

根据我的统计 meth 课程，您可以使用列表和 lapply() 以及 mapply()：

simulation <- function(N,alpha,sigma,diff,mu1)
  set.seed(12345)
  Lgroup1 <- list()
  Lgroup2 <- list()
  Lgroup1 <- lapply(1:10000, function(x) Lgroup1[[x]]<-rnorm(n=N/2, mean = mu1, sd = sigma))
  Lgroup2 <- lapply(1:10000, function(x) Lgroup2[[x]]<-rnorm(n=N/2, mean = mu1+diff, sd = sigma))
  p_values <- mapply(function(x,y) t.test(x,y)$p.value,
                     x=Lgroup1,y=Lgroup2)
  prop.table(table(p_values<alpha))[["TRUE"]]

解释：

lapply() 正在替换为第 1 组和第 2 组创建对象的循环。然后使用mapply()，我们可以获得 p 值并将其存储在向量中以供未来目标使用。

【参考方案4】：

在我的机器上，Rfast 包中的 ttest2 函数比 @BenBolker 的 sim2() 函数快一点。如果您可以忍受稍微不同的随机种子初始化，您也许可以在 linux / macos 系统上使用 parallel 包（给定多个内核和足够的内存）进一步加速该功能：

library(parallel)
library(Rfast)
#> Loading required package: Rcpp
#> Loading required package: RcppZiggurat

my_t <- function(x,y) 
  vx <- var(x)
  vy <- var(y)
  nx <- length(x)
  ny <- length(y)
  stderrx <- sqrt(vx/nx)
  stderry <- sqrt(vy/ny)
  stderr <- sqrt(stderrx^2 + stderry^2)
  df <- stderr^4/(stderrx^4/(nx - 1) + stderry^4/(ny - 1))
  tstat <- (mean(x) - mean(y))/stderr
  pval <- 2 * pt(-abs(tstat), df)
  return(pval)


sim2 <- function(N=1000,alpha=0.05,sigma=1,diff=0.1,mu1=1, nsim=1e4) 
  p_values <- numeric(nsim)  ## pre-allocate loop
  set.seed(12345)
  for(i in seq(nsim)) 
    group1 <- rnorm(n=N/2, mean = mu1, sd = sigma)
    group2 <- rnorm(n=N/2, mean = mu1+diff, sd = sigma)
    p_values[i] <- my_t(group1, group2)
  
  mean(p_values<alpha) ## replace prop.table with cheaper alternative


sim5 <- function(N=1000, alpha=0.05, sigma=1, diff=0.1, mu1=1, nsim=1e4, 
                 ncores=detectCores() - 1)
  ok <- RNGkind()
  RNGkind("L'Ecuyer-CMRG")
  set.seed(12345)
  y <- mclapply(seq(nsim), function(i)
    group1 <- rnorm(n=N/2, mean = mu1, sd = sigma)
    group2 <- rnorm(n=N/2, mean = mu1 + diff, sd = sigma)
    ttest2(group1, group2)[2]
  , mc.cores = ncores, mc.set.seed = TRUE)
  RNGkind(ok[1])
  mean(unlist(y, use.names = FALSE) < alpha)


system.time( s2 <- sim2(nsim=1e5))
#>    user  system elapsed 
#>   8.214   0.196   8.074
s2
#> [1] 0.34898
system.time( s5 <- sim5(nsim=1e5))
#>    user  system elapsed 
#>  17.196   0.573   1.712
s5
#> [1] 0.35056

^{由reprex package (v0.3.0) 于 2020 年 12 月 21 日创建}