Случайная выборка из упорядоченных данных
В моделировании нам нужны упорядоченные данные, которые представляют собой случайную выборку (с заменой или без нее) размера m из полного набора данных размера n. К сожалению, упорядочение выборочных данных оказывается узким местом в наших симуляциях. Проблема в том, что выборка повторяется R раз, что приводит к сложности выполнения O(R m log(m)). Мы стремимся снизить сложность выполнения, вызывая sort() только один раз, перед всей выборкой:
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
}
Поэтому мы задаемся вопросом, можно ли отсортировать полный набор данных только один раз, а затем напрямую получить упорядоченные выборки путем выборки из упорядоченных данных (без упорядочивания индексов, возвращаемых sample()).
Есть ли у кого-нибудь предложения, как это сделать в R (нам не обязательно нужен код R, но общего подхода тоже будет достаточно)?
Если полный набор данных упорядочен и вы выполняете выборку с использованием неупорядоченных индексов, не делает ли это выборку неупорядоченной? Как пишет @PGSA, прояснить ситуацию помог бы минимальный воспроизводимый пример.
@PGSA Я добавил пример кода, описывающий наш текущий подход и узкие места.
Если вы действительно не можете позволить себе сортировать индексы, вы могли бы, я думаю, предварительно вычислить все возможные неупорядоченные выборки (т.е. n-choose-k), что было бы очень трудоемко для большого набора, а затем случайным образом выбирать один набор на каждой итерации...
Не могли бы вы предоставить хотя бы дополнительную информацию — о каких размерах m, n и R мы здесь говорим?
n варьируется от 30 до 1000, для каждого n пробуются разные значения m, а R равно 1000. Каждая комбинация повторяется N = 10^6 раз, поскольку мы оцениваем вероятность покрытия.
Тогда предварительный расчет возможных выборок не является хорошим планом - 1000 выберите 10 (например) предлагает 2,63E + 23 возможных набора выборок.
@SamR Я отредактировал код, чтобы его можно было выполнять, а также использовал реалистичные значения для n, m и R.
Давайте продолжим обсуждение в чате.
Для выборки без замены может быть достаточно all.data[sample(rep(c(TRUE, FALSE), c(m, n-m)))].
Ответы 3
Один из возможных методов — систематическая выборка отсортированных последовательностей. Например, если из 500 требуется 100 отсортированных индексов, разделите всю выборку на 100 слоев, содержащих пять индексов, и случайным образом выберите 1 из этих пяти. Полученная выборка автоматически сортируется.
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
Сравните с алгоритмом, который сортирует индексы выборки:
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
Это происходит очень быстро (и предварительно отсортировано), но теряется элемент «настоящей» случайной выборки. Это интересный подход. Надеюсь, это обеспечит баланс, который нужен ОП.
Случайная выборка может быть простой или систематической. Оба метода дают «настоящую» случайную выборку.
Да, каждый член популяции имеет равные шансы быть выбранным в выборке, но не означает ли этот подход, что некоторые комбинации никогда не могут быть выбраны? например, ни в одной выборке не будет одновременно 1 и 2, поскольку они оба находятся в страте [1:5]? Я неправильно понял? Для меня это новая территория, поэтому я очень хочу учиться.
Да, некоторые комбинации были бы невозможны. Но каждый человек с равной вероятностью будет включен в выборку, что является основным моментом случайной выборки. Систематический отбор проб используется для экономии затрат или по другим логистическим причинам. Здесь его можно использовать вместо сортировки.
Метод @lotus — это именно то, что я искал, и он решает проблему выборки без замены. Реализация собственной процедуры выборки, как вы это сделали, несколько рискованна из-за множества ошибок при генерации случайных чисел с помощью арифметики по модулю. Таким образом, я бы не стал доверять результирующим вероятностям охвата или доверительным интервалам, полученным в результате «систематической» выборки: известно, например, что выборка складным ножом, являющаяся очень систематическим подходом, не дает репрезентативного распределения выборки.
Спасибо за ответы. Поскольку самый существенный ответ скрыт в комментарии @lotus, позвольте мне резюмировать его как ответ.
Для выборки без замены индексы могут выбираться как последовательность из n логических значений (ИСТИНА, если элемент является частью выборки) вместо последовательности из m числовых индексов. Этого можно достичь в R с помощью
selected <- sample(rep(c(TRUE, FALSE), c(m, n-m)))
sample.data <- all.data[selected]
Для выборки с заменой простой алгоритм будет следующим:
- Инициализируйте массив счетчиков c размера n нулями.
- Выберите m раз число i из {1,2,...,n} и увеличьте c[i] для каждого числа.
- В конце концов выберите каждый элемент так часто, как он учитывается в c.
Простая, хотя и довольно медленная из-за цикла for, реализация в R будет следующей:
c <- numeric(n)
is <- sample(1:n, m, replace=TRUE)
for (i in is) c[i] <- c[i] + 1
sample.data <- rep(all.data, c)
Как было предложено в комментарии @jblood94, цикл for можно обойти с помощью tabulate():
c <- tabulate(sample(n, m, replace=TRUE), n)
sample.data <- rep(all.data, c)
То же самое работает и с replace=FALSE.
Используйте tabulate, чтобы избежать цикла for: rep(x, tabulate(sample(n, m, 1), n)).
@ jblood94 Да, это действительно приводит к некоторому ускорению. Интересно, что это быстрее, чем выборка без замены из логического вектора.
Для выборки с заменой (ссылки и проверку моделирования см. ниже):
frexp <- function(m, R, x) {
y <- rexp(R)
for (r in 1:R) {
cs <- cumsum(rexp(m))
sample.data <- x[ceiling(cs*length(x)/(cs[m] + y[r]))]
}
}
Для отбора проб без замены:
fbool2 <- function(m, R, x) {
n <- length(x)
b <- logical(n)
for (r in 1:R) {
b[i <- sample(n, m)] <- TRUE
sample.data <- x[b]
b[i] <- FALSE
}
}
Другие функции для бенчмаркинга, в том числе тот, который не требует сортировки:
fNoSort <- function(m, R, x) for (r in 1:R) sample.data <- sample(x, m)
fsort <- function(m, R, x) for (r in 1:R) sample.data <- sort(sample(x, m))
fbool1 <- function(m, R, x) { # from @lotus
for (r in 1:R) sample.data <- x[sample(rep.int(!0:1, c(m, length(x) - m)))]
}
Бенчмаркинг:
microbenchmark::microbenchmark(
fNoSort = fNoSort(m, R, all.data),
fsort = fsort(m, R, all.data),
fbool1 = fbool1(m, R, all.data),
fbool2 = fbool2(m, R, all.data),
frexp = frexp(m, R, all.data)
)
#> Unit: milliseconds
#> expr min lq mean median uq max neval
#> fNoSort 8.3386 8.52895 9.746159 8.80505 9.37925 19.4995 100
#> fsort 35.6090 38.05500 40.689749 39.67485 41.73925 79.2346 100
#> fbool1 29.2213 29.64240 31.838542 30.65475 32.70145 40.3151 100
#> fbool2 10.2221 10.43100 11.616368 10.79390 12.22930 19.7826 100
#> frexp 5.0672 5.16110 6.028450 5.30465 5.81935 16.2895 100
frexp быстрее, чем sample без сортировки.
Быстрое моделирование для проверки того, что индексы, используемые в frexp, представляют собой (отсортированную) однородную случайную выборку с заменой:
(seed <- sample(.Machine$integer.max, 1))
#> [1] 904968452
set.seed(seed)
n <- 20 # sample the numbers 1 through 20
N <- 1e6 # number of samples
tabulate(sample(n, N, 1), n)
#> [1] 50176 50058 50270 50279 50208 49722 50005 49668 49856 50196 49890 49798
#> [13] 50163 49489 49918 49698 50052 50240 50333 49981
tabulate(ceiling((cs <- cumsum(rexp(N)))*n/(cs[N] + rexp(1))), n)
#> [1] 49918 49883 49957 50063 49815 50226 49792 49970 49879 49788 50119 49965
#> [13] 50417 49832 50059 50281 49930 49904 50416 49786
Использованная литература:
https://stackoverflow.com/a/63430876/9463489
Результаты производительности впечатляют. Не могли бы вы объяснить, как работает frexp() и почему генерируются случайные числа экспоненциального распределения?
Я немного подправил frexp, чтобы увеличить скорость, добавил проверку симуляции и несколько ссылок внизу.
Если бы вы могли предложить воспроизводимый пример, это очень помогло бы в выявлении и устранении узкого места.