Как сделать симуляцию без цикла?

Основываясь на моей статистике курсов метамфетамина, вы можете использовать списки и 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 и сохранить их в векторе для будущих целей.

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 против 0,0,35215 для 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

На самом деле вы можете полностью отказаться от цикла for, создав отдельные потоки group1 и group2, затем задав их размеры в виде матриц с ncol = 10000, а затем sapply для обработки t.tests:

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)

На моей машине функция ttest2 из пакета Rfast работает немного быстрее, чем функция sim2() @BenBolker. Если вы можете жить с немного другой инициализацией случайного начального числа, вы можете еще больше ускорить функцию, используя пакет parallel (с учетом нескольких ядер и достаточного объема памяти) в системе linux/macos:

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

^{Created on 2020-12-21 by the reprex package (v0.3.0)}