Подтвердить что ты не робот

R: sample() подчиняется ограничению

Я пытаюсь случайным образом пробовать 7 чисел от 0 до 7 (с заменой), но при условии ограничения, что выбранные числа складываются до 7. Так, например, вывод 0 1 1 2 3 0 0 в порядке, но выход 1 2 3 4 5 6 7 нет. Есть ли способ использовать команду sample с добавленными ограничениями?

Я намерен использовать функцию replicate() с помощью команды sample в качестве аргумента, чтобы вернуть список из N разных векторов из команды sample. То, как я в настоящее время использую команду образца (без каких-либо ограничений), мне нужно, чтобы N был очень большим, чтобы получить как можно больше возможных векторов, которые суммируются ровно до 7. Я считаю, что должен быть более простой способ сделать это!

Вот мой код для этой части:

x <- replicate(100000, sample(0:7, 7, replace=T))    

В идеале я хочу, чтобы 10 000 или 100 000 векторов в x суммировались до 7, но для этого понадобилось бы огромное значение N. Спасибо за любую помощь.

4b9b3361

Ответ 1

Чтобы убедиться, что вы выборочно выполняете выборку, вы можете просто сгенерировать все перестановки и ограничить те, которые суммируются с 7:

library(gtools)
perms <- permutations(8, 7, 0:7, repeats.allowed=T)
perms7 <- perms[rowSums(perms) == 7,]

Из nrow(perms7), мы видим, что существует только 1716 возможных перестановок, суммируемых в 7. Теперь вы можете равномерно выбирать из перестановок:

set.seed(144)
my.perms <- perms7[sample(nrow(perms7), 100000, replace=T),]
head(my.perms)
#      [,1] [,2] [,3] [,4] [,5] [,6] [,7]
# [1,]    0    0    0    2    5    0    0
# [2,]    1    3    0    1    2    0    0
# [3,]    1    4    1    1    0    0    0
# [4,]    1    0    0    3    0    3    0
# [5,]    0    2    0    0    0    5    0
# [6,]    1    1    2    0    0    2    1

Преимущество такого подхода состоит в том, что легко видеть, что мы выборочно выбираем случайным образом. Кроме того, он довольно быстро - здание perms7 заняло 0,3 секунды на моем компьютере и построило 1 миллион строк my.perms заняло 0,04 секунды. Если вам нужно нарисовать много векторов, это будет немного быстрее, чем рекурсивный подход, потому что вы просто используете индексирование матриц в perms7 вместо генерации каждого вектора отдельно.

Здесь распределение отсчетов чисел в выборке:

#      0      1      2      3      4      5      6      7 
# 323347 188162 102812  51344  22811   8629   2472    423 

Ответ 2

Начните со всех нулей, добавьте их в любой элемент, выполните 7 раз:

sumTo = function(){
    v = rep(0,7)
    for(i in 1:7){
        addTo=sample(7)[1]
        v[addTo]=v[addTo]+1
    }
    v
}

Или, что то же самое, просто выберите, какой из 7 элементов вы собираетесь увеличивать в одном примере длины 7, а затем табулируйте их, убедившись, что вы вставляете в таблицу до 7:

sumTo = function(){tabulate(sample(7, 7, replace = TRUE), 7)}


> sumTo()
[1] 2 1 0 0 4 0 0
> sumTo()
[1] 1 3 1 0 1 0 1
> sumTo()
[1] 1 1 0 2 1 0 2

Я не знаю, приведет ли это к единому образцу из всех возможных комбинаций...

Распределение отдельных элементов более 100 000 повторений:

> X = replicate(100000,sumTo())
> table(X)
X
     0      1      2      3      4      5      6 
237709 277926 138810  38465   6427    627     36 

Не поражал 0,0,0,0,0,7 того времени!

Ответ 3

Этот рекурсивный алгоритм выведет распределение с большей вероятностью для больших чисел, чем другие решения. Идея состоит в том, чтобы сбросить случайное число y в 0:7 в любом из семи доступных слотов, затем повторить со случайным числом в 0:(7-y) и т.д.:

sample.sum <- function(x = 0:7, n = 7L, s = 7L) {
   if (n == 1) return(s)
   x <- x[x <= s]
   y <- sample(x, 1)
   sample(c(y, Recall(x, n - 1L, s - y)))
}

set.seed(123L)
sample.sum()
# [1] 0 4 0 2 0 0 1

Рисунок 100 000 векторов занял 11 секунд на моей машине, и вот распределение, которое я получаю:

#      0      1      2      3      4      5      6      7 
# 441607  98359  50587  33364  25055  20257  16527  14244 

Ответ 4

Может быть более простой и/или более элегантный способ, но здесь метод грубой силы, использующий функцию LSPM:::.nPri. Ссылка включает определение для версии алгоритма только для R. Для тех, кого это интересует.

#install.packages("LSPM", repos="http://r-forge.r-project.org")
library(LSPM)
# generate all possible permutations, since there are only ~2.1e6 of them
# (this takes < 40s on my 2.2Ghz laptop)
x <- lapply(seq_len(8^7), nPri, n=8, r=7, replace=TRUE)
# set each permutation that doesn't sum to 7 to NULL
y <- lapply(x, function(p) if(sum(p-1) != 7) NULL else p-1)
# subset all non-NULL permutations
z <- y[which(!sapply(y, is.null))]

Теперь вы можете выбрать из z и быть уверены, что получаете перестановку, суммируемую до 7.

Ответ 5

Я задаю этот вопрос интригующим и дал ему дополнительную мысль. Другой (более общий) подход к (приближенному) образцу равномерно из всех возможных решений без генерации и хранения всех перестановок (что явно не представляется возможным в случае с гораздо более чем 7 числами), в R на sample() может быть простая реализация MCMC:

S <- c(0, 1, 1, 2, 3, 0, 0) #initial solution
N <- 100 #number of dependent samples (or burn in period)
series <- numeric(N)
for(i in 1:N){
    b <- sample(1:length(S), 2, replace=FALSE) #pick 2 elements at random
    opt <- sum(S[-b]) #sum of complementary elements
    a <- sample(0:(7-opt), 1) #sample a substistute
    S[b[1]] <- a #change elements
    S[b[2]] <- 7 - opt - a 
}
S #new sample 

Это, конечно, очень быстро для нескольких образцов. "Распределение":

#"distribution" N=100.000:      0      1      2      3      4      5      6      7
#                            321729 189647 103206  52129  22287   8038   2532    432

Конечно, в этом случае, когда на самом деле можно найти и сохранить все комбинации, и если вы хотите получить огромную выборку из всех возможных результатов, просто используйте partitions::compositions(7, 7), также предложенный Джошем О'Брайеном в комментариях, чтобы избежать вычисления всех перестановок, когда требуется лишь малая доля:

perms7 <- partitions::compositions(7, 7)

>tabulate(perms7[, sample(ncol(perms7), 100000, TRUE)]+1, 8)
#"distribution" N=100.000:      0      1      2      3      4      5      6      7
#                            323075 188787 102328  51511  22754   8697   2413    435