Параллельное моделирование Монте-Карло в R с использованием снегопада

Я пытаюсь сравнить до тысяч предполагаемых бета-дистрибутивов. Каждое бета-распределение характеризуется двумя параметрами формы альфа и бета. Сейчас я рисую по 100 000 образцов каждого дистрибутива. В качестве окончательного результата я хочу получить порядок распределений с наивысшей вероятностью в каждом розыгрыше выборки. Мой первый подход заключался в использовании lapply для создания матрицы числовых значений N * NDRAWS, которая потребляла слишком много памяти, поскольку N превышала 10 000. (10,000 * 100,000 * 8 байтов)

Поэтому я решил использовать последовательный подход к упорядочиванию каждого розыгрыша, затем суммировать порядок всех розыгрышей и получить окончательный порядок, как показано в примере ниже:

set.seed(12345)
N=100
NDRAWS=100000
df <- data.frame(alpha=sample(1:20, N, replace=T), beta=sample(1:200, N, replace=T))  

vec    <- vector(mode = "integer", length = N )

for(i in 1:NDRAWS){
  # order probabilities after a single draw for every theta
  pos <- order(rbeta(N, shape1=df$alpha, shape2=df$beta) )

  # sum up winning positions for every theta
  vec[pos] <- vec[pos] + 1:N
}

# order thetas
ord <- order(-vec)

df[ord,]

Это потребляет только N * 4 байта памяти, так как нет гигантской матрицы, а есть единственный вектор длины N. Мой вопрос теперь в том, как ускорить эту операцию с помощью снегопада (или любого другого многоядерного пакета), воспользовавшись моим преимуществом. 4 ядра процессора вместо одного ядра ???

# parallelize using snowfall pckg
library(snowfall)
sfInit( parallel=TRUE, cpus=4, type="SOCK")
sfLapply( 1:NDRAWS, function(x) ?????? )
sfStop()

Любая помощь приветствуется!