Efficiently Sampling Lottery Numbers in R

user_12576489 0
Viewed 710

I want to write a function that samples n lottery tickets with 6 numbers each from 1 to 45 without replacement. However, I need to do this efficiently, meaning no loops or loop-like functions. (I guess Rcpp would also work but I would prefer a vectorized solution in base R)

Solving without the restriction:

lottery_inef <- function(n){
  
 t(replicate(n,
          sample(1:45, 6)))
}

So here I get a matrix where each row corresponds to one lottery ticket. Now if I want to simulate millions of lottery tickets this gets pretty slow and hence I am interested in a vectorized solution.

My Idea here would be:

lottery_ef <- function(n){
  
  m <- matrix(sample(1:45, n*6, replace = TRUE), ncol = 6)
  
  # somehow subset the matrix without a loop to remove all the 
  # rows that have non-unique values as in the lottery we can only draw each number once
}

For the efficient version, I am sort of lost at the point of subsetting without a loop or apply(). If anyone can solve this subsetting problem or point me in a completely different direction that will lead me to a solution I would be very grateful.

sampling r function
3 Answers

replicate doesn't actually do all that well at this scale. With Just In Time compilation (used in R for several years now) a for loop can be faster, especially when we can precisely pre-allocate the data structure. We can also avoid the t():

lottery_inef <- function(n){
 t(replicate(n,
          sample(1:45, 6)))
}

lottery_preall <- function(n){
  m = matrix(NA_integer_, nrow = n, ncol = 6)
  for(i in 1:n) {
    m[i, ] = sample.int(45L, size = 6)
  }
  m
}

nn = 1e6
microbenchmark::microbenchmark(
  lottery_inef(nn), 
  lottery_preall(nn),
  times = 2
)
# Unit: seconds
#                expr      min       lq     mean   median       uq      max neval
#    lottery_inef(nn) 9.400862 9.400862 9.571756 9.571756 9.742649 9.742649     2
#  lottery_preall(nn) 4.948216 4.948216 5.454482 5.454482 5.960749 5.960749     2

replicate accumulates results in a list and then needs to check the dimensions of each one before deciding it can be simplified to a matrix, and has to do that conversion. All of that overhead is skipped with a pre-allocated integer matrix for a roughly 2x speed-up.

We could also compare to, say vapply (a quick test shows vapply just a little slower than the loop), but I think to get much more speed out of this you'd need to run in parallel - which would be a good option here and could probably get you a speed-up almost equal to the number of cores you use.

sample.int is pretty much just a call to C code, so going to Rcpp probably won't do that much better - I think parallelization is your best bet for more speed.

user_9030610

Since generating all of the combinations for a set this size only takes a few seconds, it might be worthwhile to do that and then subset it for the 'lottery tickets'. Below I've used sample() to generate 1 million row indices (both with and without replacement), and bracketted subsetting on the entire set to generate possible tickets.

If you need to do this often enough, or at different times it might be worthwhile to save the full combination set rather than re-generate it each time. Nearly all of the processing is generating the full set of combinations. Selecting the 'tickets' after that is quick.

Timings show it took ~6 seconds to create all the combinations, ~.2 seconds for 1 million indices, and ~.1 sec for bracketted subsetting of 1 million rows.

set.seed(2)

tictoc::tic() #included for timing

# All possible lotto combinations as matrix, 1 per row
lotto_all <- t(combn(1:45, 6))

tictoc::toc() #included for timing
#> 5.899 sec elapsed

# A look at the data:
head(lotto_all)
#>      [,1] [,2] [,3] [,4] [,5] [,6]
#> [1,]    1    2    3    4    5    6
#> [2,]    1    2    3    4    5    7
#> [3,]    1    2    3    4    5    8
#> [4,]    1    2    3    4    5    9
#> [5,]    1    2    3    4    5   10
#> [6,]    1    2    3    4    5   11

# Getting index (row) numbers for our 'tickts' with & without replacement
tictoc::tic()
sample_indices_no_replacement <- sample(1:nrow(lotto_all), size = 1e6, replace = F)
tictoc::toc()
#> 0.178 sec elapsed

sample_indices_w_replacement <- sample(1:nrow(lotto_all), size = 1e6, replace = T)

# The number combinations of our 'tickets'
tictoc::tic()
sample_tickets_no_rep <- lotto_all[sample_indices_no_replacement,]
tictoc::toc()
#> 0.097 sec elapsed

sample_tickets_rep <- lotto_all[sample_indices_w_replacement,]

# A look at the sample tickets:
head(sample_tickets_no_rep)
#>      [,1] [,2] [,3] [,4] [,5] [,6]
#> [1,]    8   12   14   31   34   44
#> [2,]    6   10   16   26   32   36
#> [3,]    3    4   10   15   41   43
#> [4,]    2    3    5   17   33   36
#> [5,]    7   17   24   25   35   40
#> [6,]   32   33   34   36   39   43

# See that there are some duplicates using replacement = T
length(unique(sample_indices_no_replacement))
#> [1] 1000000
length(unique(sample_indices_w_replacement))
#> [1] 941309

Created on 2020-10-27 by the reprex package (v0.3.0)

user_75473270

Since efficiency is the main focus, there are a couple of packages, arrangements and RcppAlgos*, that have this specific task in mind.

Before we get started, we first illustrate that when we are using sample, we have no control over the uniqueness of the outcome. Each draw is from a uniform distribution, thus it is possible that we could repeatedly draw the same permutation multiple times. Using the functions by @Gregor, we have:

set.seed(42)
system.time(a <- lottery_inef(1e6))
 user  system elapsed 
7.640   0.345   7.984

sum(duplicated(a))
[1] 86

set.seed(42)
system.time(b <- lottery_preall(1e6))
 user  system elapsed 
3.673   0.256   3.929

sum(duplicated(b))
[1] 86

While using the package arrangements is faster, we still see the same behavior:

set.seed(42)
system.time(arng <- arrangements::permutations(45, 6, nsample = 1e6))
  user  system elapsed 
 0.761   0.021   0.785 

sum(duplicated(arng))
[1] 108

Now, using the package RcppAlgos we are guaranteed that the result is unique if the number of results requested is less than the total number of results (over 5 billion in our case):

RcppAlgos::permuteCount(45, 6)
[1] 5864443200

system.time(algosSer <- RcppAlgos::permuteSample(45, 6,
                                                 n = 1e6,
                                                 seed = 42))
 user  system elapsed 
0.560   0.001   0.561

sum(duplicated(algosSer))
[1] 0

Additionally, we can utilize multiple threads via the nThreads argument for an even greater speed up.

system.time(algosPar <- RcppAlgos::permuteSample(45, 6,
                                                 n = 1e6,
                                                 seed = 42,
                                                 nThreads = 4))
 user  system elapsed 
0.574   0.001   0.280

## Results are the same as the serial version
identical(algosPar, algosSer)
[1] TRUE

* I am the author of RcppAlgos

user_44085380
Related