Random sampling from ordered data

Question

In a simulation, we need ordered data which is a random sample (with or without replacement) of size m from a full data set of size n. Unfortunately, ordering the sampled data turns out to be a performance bottleneck in our simulations. The problem is that the sampling is repeated R times, which results in a runtime complexity of O(R m log(m)). We are seeking to reduce the runtime complexity by calling sort() only once, before all sampling:

n <- 500; m <- 100; R <- 1000
all.data <- runif(n)
all.data <- sort(all.data)  # the full data set is already sorted
for (r in 1:R) {
  indices <- sample(1:n, m)
  sample.data <- sort(all.data[indices])  # this call to sort should be avoided
}

We therefore wonder whether it is possible to sort the full data set only once and then subsequently to directly obtain ordered samples by sampling from the ordered data (without ordering the indices returned by sample()).

Does anyone have a suggestion how to do this in R (we do not necessarily need R code, but the general approach would also be sufficient)?

Answer 1

One possible method is to do systematic sampling on sorted sequences. For example, if 100 sorted indices are desired from 500, then partition the entire sample into 100 strata containing five indices and randomly sample 1 of these five. The resulting sample is automatically sorted.

pseudo.sort <- function() {
  size <- n/m
  indices <- sapply(1:m, function(i) sample(1:size + (size * (i-1)), 1))
  replicate(R, \() all.data[indices])
}

R <- 10000

pseudo.sort()
  [1]   4   8  13  17  21  29  35  38  44  46  55  56  61  69  72  77  82  89  93  98
 [21] 102 109 115 120 123 127 132 136 143 149 154 156 162 167 171 180 185 186 194 198
 [41] 202 208 214 216 221 226 233 240 243 246 255 257 264 266 271 276 284 290 294 300
 [61] 301 307 314 316 325 326 333 340 341 347 353 357 365 367 374 376 384 388 395 397
 [81] 404 406 414 418 422 428 432 440 445 447 454 456 464 468 471 477 485 489 492 496

Compare with an algorithm which sorts the sample indices:

sorted.sample <- \() {
  for (r in 1:R) {
    indices <- sort(sample(1:n, m)) # the sort has been shifted
    sample.data <- all.data[indices]  
  }
}

lotus.sort <- \() replicate(R, \() all.data[sample(rep(c(TRUE, FALSE), c(m, n-m)))])

> microbenchmark::microbenchmark(sorted.sample(), pseudo.sort(), lotus.sort())
Unit: microseconds
            expr     min      lq      mean   median       uq      max neval
 sorted.sample() 44166.7 45439.5 50490.146 45936.05 47094.55 363412.8   100
   pseudo.sort()  1166.3  1201.5  1232.323  1225.90  1258.80   1400.6   100
    lotus.sort()   535.3   555.00   612.407   569.55   586.65  4126.2   100