Более эффективный способ подсчета доли случайных чисел ниже порога

Question

Работая над некоторыми вероятностными упражнениями, мне нужно изобразить долю случаев, когда монета выпадает в голову (например, недобросовестная монета с р = 0,3) - против количества подбрасываний монеты.

Здесьмой вдохновленный Python-код R (он компилируется и выполняется), который работает очень медленно.Есть ли способ сделать его более идиоматическим R-кодом?

Очень ценится

flip_experiment <- function(prob_heads, n_flips) {
  proportion_heads <- c()
  for (i in 1:n_flips) {
    count = 0
    for (j in 1:i){
      if(runif(1, 0, 1) <= prob_heads){
        # We flipped a head!
        count <- count + 1
      }
    }
    prop_heads = count / i
    proportion_heads <- append(proportion_heads, prop_heads)

  }
  plot(1:n_flips, proportion_heads)
  return 
}

flip_experiment(0.3, 1000);
flip_experiment(0.03, 1000);

Вот один из графиков:

Answer 1

Решение на основе

a tidyverse также очень эффективно:

library(tidyverse)
flip.experiment <- function(prob.heads, n.flips){
  dat = tibble(exp.no = 1:n.flips, flip.no = 1:n.flips)
  dat = dat %>% 
    group_by(exp.no) %>% 
    expand(flip.no = 1:exp.no) %>%
    mutate(head = if_else(runif(exp.no, 0, 1) <= prob.heads, TRUE, FALSE)) %>%
    summarise(proportion.heads = sum(head)/n())
  plot <- plot(dat$exp.no, dat$proportion.heads)
  return(plot)
}
flip.experiment(0.3, 1000)

Answer 2

R прекрасно, потому что в большинстве случаев вам не нужны петли:

set.seed(1)
x <- 1:1000
y <- sapply(x, function(j) sum(rbinom(j, size = 1, prob = .3))/j)
plot(x, y)

Пояснение

Бросок монеты следует биномиальному распределению.В R вы можете получить случайные результаты для ряда распределений (см. ?distribution), одно из которых является биномиальным с rbinom (здесь r означает случайное).Бит sapply используется для упрощенного применения, векторизации.За sum и делением следите наверняка;)

Answer 3

Вот функция, которую я предложил, используя функцию rbinom.

flip_experiment2 <- function(prob_heads, n_flips){
  proportion_heads <- rbinom(n_flips, 1:n_flips, prob_heads)/(1:n_flips)
  return(plot(x = 1:n_flips, y = proportion_heads))
}

Здесь я провел анализ производительности.Использование rbinom намного быстрее.

library(microbenchmark)

microbenchmark(m1 = flip_experiment(0.3, 1000), 
               m2 = flip_experiment2(0.3, 1000))

Unit: milliseconds
 expr       min        lq      mean    median        uq       max neval cld
   m1 765.22263 859.28831 923.04026 900.44548 970.77151 1259.4624   100   b
   m2  27.88089  29.93146  33.50223  31.55485  33.50544  146.7657   100  a 

Answer 4

По крайней мере j -loop можно векторизовать (используя runif(i, 0, 1) вместо вызова runif(1, 0, 1) для i раз).Кроме того, не увеличивайте результирующий вектор proportion_heads.Его размер известен, так что вы можете предварительно выделить его и заполнить.

flip_experiment <- function(prob_heads, n_flips) {
  proportion_heads <- numeric(n_flips)
  for (i in 1:n_flips) {
    count <- sum(runif(i, 0, 1) <= prob_heads)
    proportion_heads[i] <- count / i  
    }
  plot(1:n_flips, proportion_heads)
  }

flip_experiment(0.3, 1000);
flip_experiment(0.03, 1000);

Вторая мысль:

flip_experiment <- function(prob_heads, n_flips) {
  proportion_heads <- cumsum(runif(n_flips, 0, 1) <= prob_heads) / seq_len(n_flips)
  plot(1:n_flips, proportion_heads)
  }