выборка (.) в сравнении с рассчитанными вероятностями

Question

Где моя ошибка в следующем коде / аргументации.

Если у меня есть n предметов, и я хочу посчитать количество комбинаций (т.е. порядок не имеет значения) розыгрышей размером k с заменой , я могу использовать биномиальный коэффициент:

choose(n + k - 1, k)

Когда я хочу посчитать количество возможностей, в которых элемент i не является частью розыгрыша, я бы, таким образом, использовал

choose(n + k - 1 - 1, k)

(В принципе, у меня есть на один вариант меньше, чтобы рисовать).

например. с n=3, k= 3 у меня есть следующие 10 (choose(3 + 3 - 1, 3)) розыгрыши:

# AAA, AAB, ABB, BBB, BBC, BCC, CCC, CCA, CAA, ABC 

Те, где элемент 'A' (скажем) не появляется, это следующие 4 (choose(3 + 3 - 1 - 1, 3)):

# BBB, BBC, BCC, CCC

Пока все хорошо. Таким образом, я могу вычислить вероятность того, что у меня есть ничья с n = k, где элемент i не появляется:

freeFromCounts <- function(n) choose(2 * n - 2, n)
totalCounts    <- function(n) choose(2 * n - 1, n)
ratio          <- function(n) freeFromCounts(n) / totalCounts(n) ## (n - 1) / (2 * n - 1)

Итак, вот моя проблема, если я имитирую розыгрыши (с sample(.)) и повторяю это n.rep раз, я ожидаю увидеть около n.rep * ratio(n) розыгрышей, где элемент 1 (скажем) отсутствует. Это, однако, не тот случай. Где моя ошибка?

sim <- function(n, n.rep = 10000, x0 = 1) {
    sum(replicate(n.rep, {
       s <- sample(n, n, TRUE)
       all(s != x0)
    })) / n.rep
}

set.seed(12)
sim(10, 1e6)  # [1] 0.348278
ratio(10)     # [1] 0.4736842

Answer 1

Коллега нашел решение. sample сэмплы не неупорядоченные , но заказанные , т. Е. Разница между [1, 2, 3] и [3, 2, 1]. С этим знанием цифры одинаковы: 9 ^ 10 / 10 ^ 10 = 0.3486784.

Итак, чтобы выполнить симуляцию для исходной задачи, вам нужно придумать вариант sample, который делает неупорядоченную выборку . Идея заключается в известном выводе формулы для неупорядоченной выборки с заменой (см., Например, http://mathworld.wolfram.com/Multichoose.html)

sampleUnordered <- function(n, k) {
  aS <- rep(0, n + k -1)
  aS[!seq_along(aS) %in% sample(n + k - 1, k)] <- seq.int(n-1)
  rl <- rle(aS)
  zeros <- which(rl$values == 0)
  good <- rl$lengths[zeros]
  nr <- c(rl$values, n)[zeros + 1]
  rep(nr, good)
}

sim <- function(n, n.rep = 10000, x0 = 1) {
    sum(replicate(n.rep, {
       s <- sampleUnordered(n, n)
       all(s != x0)
    })) / n.rep
}

set.seed(12)
sim(10, 1e6) # [1] 0.473234