Каков самый быстрый способ R для генерации случайных переменных распределения Бейтса?

Получить случайные числа распределения Бейтса несложно. Мне нужно 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() ).

На моем старом ноутбуке 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)

Gregor Thomas 23.05.2023 06:25
Стоит ли изучать PHP в 2023-2024 годах?
Стоит ли изучать PHP в 2023-2024 годах?
Привет всем, сегодня я хочу высказать свои соображения по поводу вопроса, который я уже много раз получал в своем сообществе: "Стоит ли изучать PHP в...
Поведение ключевого слова "this" в стрелочной функции в сравнении с нормальной функцией
Поведение ключевого слова "this" в стрелочной функции в сравнении с нормальной функцией
В JavaScript одним из самых запутанных понятий является поведение ключевого слова "this" в стрелочной и обычной функциях.
Приемы CSS-макетирования - floats и Flexbox
Приемы CSS-макетирования - floats и Flexbox
Здравствуйте, друзья-студенты! Готовы совершенствовать свои навыки веб-дизайна? Сегодня в нашем путешествии мы рассмотрим приемы CSS-верстки - в...
Тестирование функциональных ngrx-эффектов в Angular 16 с помощью Jest
В системе управления состояниями ngrx, совместимой с Angular 16, появились функциональные эффекты. Это здорово и делает код определенно легче для...
Концепция локализации и ее применение в приложениях React ⚡️
Концепция локализации и ее применение в приложениях React ⚡️
Локализация - это процесс адаптации приложения к различным языкам и культурным требованиям. Это позволяет пользователям получить опыт, соответствующий...
Пользовательский скаляр GraphQL
Пользовательский скаляр GraphQL
Листовые узлы системы типов GraphQL называются скалярами. Достигнув скалярного типа, невозможно спуститься дальше по иерархии типов. Скалярный тип...
1
1
56
4
Перейти к ответу Данный вопрос помечен как решенный

Ответы 4

Ответ принят как подходящий

В 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) } как вы предложили

bprijadi 24.05.2023 08:57

извините, мне пришлось явно указать библиотеку, такую ​​​​как { times(10000) %dopar% (sum(dqrng::dqrunif (1e6, 1, 5)) / 1e6) }, чтобы работать

bprijadi 24.05.2023 09:09

Здесь вы не моделируете распределение Бейтса. Распределение Бейтса с целым параметром 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)

Другие вопросы по теме