Получить случайные числа распределения Бейтса несложно. Мне нужно 1 миллион средних значений для запуска 10 000 раз:
bates1 = replicate(10000, mean(runif (1e+06,1,5)))
summary(bates1)
Я ждал целую вечность, чтобы завершить его расчеты. Я пробовал циклы тоже безрезультатно (бесконечно медленно).
Любой выход из этого?
Я попробовал цикл for,
set.seed(999)
for (i in 1:10000) {
x <- randomLHS(1e+6,1)
x <- 1 + 4*x
y[i] <- mean(x)
}
summary(y)
И перед кодом выделение места для x и y (используя length() ).
В R существует множество способов параллельных вычислений. Вы можете посмотреть:
В качестве примера использую библиотеку doParallel
на моей рабочей машине (скромный Surface Book 2):
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 минут (не совсем «вечно», но достаточно долго).
Некоторые из этих других ответов также могут помочь.
Без распараллеливания мы можем работать быстрее, находя высокопроизводительные версии необходимых функций. Пакет dqrng имеет версию runif
, которая примерно в 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 минутами для исходного кода на моем ноутбуке).
Спасибо за помощь, Грегор, используя точный код @neilfws (7 ядер), я получил 88 секунд system.time( { times(10000) %dopar% mean(runif (1e+06,1,5)) } )
, но не смог запустить { times(10000) %dopar% (sum(dqrunif (1e6, 1, 5)) / 1e6) }
как вы предложили
извините, мне пришлось явно указать библиотеку, такую как { times(10000) %dopar% (sum(dqrng::dqrunif (1e6, 1, 5)) / 1e6) }
, чтобы работать
Здесь вы не моделируете распределение Бейтса. Распределение Бейтса с целым параметром n является средним значением n равномерных случайных величин на (0,1). Быстрый способ получить это:
n <- 6
nsims <- 100000
usims <- matrix(runif (n*nsims), nrow = n, ncol = nsims)
bates_sims <- colMeans(usims)
Распределение Бейтса можно рассматривать как полиномиальное приближение к нормальному распределению (когда 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)
На моем старом ноутбуке
bench::mark(mean(runif (1e6, 1, 5)))
показывает 25,5 итераций в секунду, поэтому 10 тыс. итераций должно занимать около 6,5 минут. Я запустил цикл (следующий), который занял 6,8 минуты, поэтому было около 20 секунд дополнительных накладных расходов. Все это кажется прекрасным - 25,5 итераций в секунду не кажутся "бесконечно медленными". Вот мой код цикла:start = Sys.time(); n = 1e4; bates = numeric(n); for(i in 1:n) { bates[i] = mean(runif (1e6, 1, 5)) }; (run = Sys.time() - start)