How to efficiently apply the rbinom function to each row in a data frame?

Question

Given a data table with the count and rates of change for different variables, how can I sample from the count of each variable given the rate? For example given the following data table, I can loop through and use the sample or rbinorm function to get the desired output. However, the dataset I am trying to implement this on is very large. Is there a method to improve performance?

library(data.table)
set.seed(1)

dt <- data.table(
count = sample(10000:20000, 100),
rate = sample(1:20, 100, replace = T) / 1000
)

system.time(
for (i in 1:nrow(dt)){
  dt$sample_n[i] <- sum(sample(1:0, 
                           dt$count[i], 
                           prob = c(dt$rate[i], 1-dt$rate[i]), 
                           replace = T))
}
)

system.time(
for (i in 1:nrow(dt)){
  dt$sample_n2[i] <- rbinom(size = dt$count[i], n = 1, prob = dt$rate[i])
}
)

Answer 1

Assign by reference with := with just rbinom() (no loops).

Setup

library(data.table)
options(datatable.print.class = TRUE)

sample_size <- 5e4

dt <- data.table(
  count = sample(seq(10000, 10000 + sample_size), size = sample_size),
  rate = sample(1:20, size = sample_size, replace = TRUE) / 1000
)

Solution

dt[, sample_n := rbinom(n = .N, size = dt$count, prob = rate)]
dt
#>        count  rate sample_n
#>        <int> <num>    <int>
#>     1: 26100 0.016      431
#>     2: 15145 0.008      114
#>     3: 24952 0.001       23
#>     4: 31437 0.020      621
#>     5: 58358 0.008      468
#>    ---                     
#> 49996: 30517 0.002       56
#> 49997: 59047 0.009      500
#> 49998: 48737 0.018      896
#> 49999: 29686 0.005      152
#> 50000: 52429 0.011      580

Shootout

results <- list()

set.seed(1)
dt <- data.table(
  count = sample(seq(10000, 10000 + sample_size), size = sample_size),
  rate = sample(1:20, size = sample_size, replace = TRUE) / 1000
)

Original Timing

results$original1_no_modify <- system.time( # not modifying `dt`
for (i in 1:nrow(dt)) {
  sum(
    sample(1:0, dt$count[i], prob = c(dt$rate[i], 1L - dt$rate[i]), replace = TRUE)
  )
}
)

set.seed(1)
results$original1_modify <- system.time( # modifying `dt`
for (i in 1:nrow(dt)) {
  dt$sample_n[i] <- sum(
    sample(1:0, dt$count[i], prob = c(dt$rate[i], 1L - dt$rate[i]), replace = TRUE)
  )
}
)


results$original2_no_modify <- system.time( # not modifying `dt`
for (i in 1:nrow(dt)){
  rbinom(size = dt$count[i], n = 1L, prob = dt$rate[i])
}
)

set.seed(1)
results$original2_modify <- system.time( # modifying `dt`
for (i in 1:nrow(dt)){
  dt$sample_n2[i] <- rbinom(size = dt$count[i], n = 1L, prob = dt$rate[i])
}
)

:= + mapply() + rbinom() (faster, but still R-level iteration)

results$mapply_no_modify <- system.time( # not modifying `dt`
mapply(
  function(.count, .rate) rbinom(size = .count, n = 1L, prob = .rate),
  dt$count, dt$rate 
)
)

set.seed(1)
results$mapply_modify <- system.time( # modifying `dt`
dt[, sample_n3 := mapply(
  function(.count, .rate) rbinom(size = .count, n = 1L, prob = .rate),
  count, rate 
)]
)

Solution

results$solution_no_modify <- system.time( # not modifing `dt`
rbinom(n = nrow(dt), size = dt$count, prob = dt$rate)
)

set.seed(1)
results$solution_modify <- system.time(
dt[, sample_n4 := rbinom(n = .N, size = dt$count, prob = rate)]
)

Final Data Frame

dt[]
#>        count  rate sample_n sample_n2 sample_n3 sample_n4
#>        <int> <num>    <int>     <int>     <int>     <int>
#>     1: 34387 0.009      295       310       310       310
#>     2: 53306 0.019     1076      1004      1004      1004
#>     3: 14049 0.019      268       247       247       247
#>     4: 21570 0.002       45        55        55        55
#>     5: 35172 0.009      313       346       346       346
#>    ---                                                   
#> 49996: 37432 0.020      724       722       722       722
#> 49997: 14985 0.006       82        76        76        76
#> 49998: 16007 0.007      107       106       106       106
#> 49999: 49298 0.003      145       140       140       140
#> 50000: 41427 0.001       49        40        40        40

Sanity Check

stopifnot(
  identical(dt$sample_n2, dt$sample_n3) &&
  identical(dt$sample_n3, dt$sample_n4)
)

Results

results
#> $original1_no_modify
#>    user  system elapsed 
#>  18.713   0.568  19.288 
#> 
#> $original1_modify
#>    user  system elapsed 
#>  28.217   0.020  28.237 
#> 
#> $original2_no_modify
#>    user  system elapsed 
#>   0.155   0.000   0.155 
#> 
#> $original2_modify
#>    user  system elapsed 
#>   9.085   0.152   9.237 
#> 
#> $mapply_no_modify
#>    user  system elapsed 
#>   0.139   0.000   0.139 
#> 
#> $mapply_modify
#>    user  system elapsed 
#>   0.132   0.000   0.131 
#> 
#> $solution_no_modify
#>    user  system elapsed 
#>   0.004   0.000   0.004 
#> 
#> $solution_modify
#>    user  system elapsed 
#>   0.004   0.000   0.004
rbindlist(lapply(results, as.list), idcol = "approach")
#>               approach user.self sys.self elapsed user.child sys.child
#>                 <char>     <num>    <num>   <num>      <num>     <num>
#> 1: original1_no_modify    18.713    0.568  19.288          0         0
#> 2:    original1_modify    28.217    0.020  28.237          0         0
#> 3: original2_no_modify     0.155    0.000   0.155          0         0
#> 4:    original2_modify     9.085    0.152   9.237          0         0
#> 5:    mapply_no_modify     0.139    0.000   0.139          0         0
#> 6:       mapply_modify     0.132    0.000   0.131          0         0
#> 7:  solution_no_modify     0.004    0.000   0.004          0         0
#> 8:     solution_modify     0.004    0.000   0.004          0         0

Answer 2

All sampling functions are usually vectorized, meaning you can directly do:

dt$sample_n2 <- rbinom(size = dt$count, n = nrow(dt), prob = dt$rate)