Это код, который я сейчас использую:
n <- 7
N <- 52
r <- 13
reps <- 1000000
deck <- rep(c('h','d','c','s'), each = r)
diamonds <- rep(NA, length.out = reps)
pos <- sample(x = 1:52, size = 7, replace = FALSE)
for(i in 1:reps) {
hand <- sample(x = deck, replace = FALSE)[pos]
diamonds[i] <- sum(ifelse(hand == 'd', 1, 0))
}
barplot(table(diamonds), col = 'red', xlab = '# of diamonds',
ylab = paste('frequency out of',reps,'trials'),
main = paste('Positions:',pos[1],pos[2],pos[3],pos[4],
pos[5],pos[6],pos[7]))
Что мне действительно нравится, так это возможность присвоить заголовок барплоту примерно так:
barplot(..., main = paste('Positions:',pos))
и иметь заголовок «Позиции: p1 p2 p3 p4 p5 p6 p7», где p1, p2, ..., p7 - элементы поз.
Для всех, кто интересуется, этот код случайным образом выбирает 7 позиций из 52, а затем подсчитывает количество алмазов ('d') в этих 7 позициях после каждой перетасовки колоды на 1000000 перетасовок. Затем строится эмпирическое распределение количества алмазов в этих 7 картах.