Я пытаюсь обернуть голову вокруг этой проблемы.
Плашка бросается 100 раз. Какова вероятность того, что лицо не появится более 20 раз? Моей первой мыслью было использование биномиального распределения P (x) = 1 - 6 смс (100, 1/6, 20), но это, очевидно, неправильно, поскольку мы считаем несколько случаев более одного раза. Моя вторая идея состоит в том, чтобы перечислить все возможные броски x1 + x2 + x3 + x4 + x5 + x6 = 100, чтобы xi <= 20, и суммировать многочлены, но это кажется слишком сложным в вычислительном отношении. Приблизительные решения также будут работать для меня.

Это обобщение известной проблемы дня рождения : если человек, у которых есть случайные, равномерно распределенные «дни рождения» среди набора возможностей, какова вероятность того, что ни один день рождения не разделен более чем человек?д = 6 м = $n=100$$d=6$$m=20$

Точный расчет дает ответ (с двойной точностью). Я набросаю теорию и приведу код для общих Асимптотическая синхронизация кода - что делает его пригодным для очень большого числа дней рождения и обеспечивает разумную производительность, пока не исчисляется тысячами. В этот момент приближение Пуассона, обсуждавшееся при распространении парадокса дня рождения на более чем 2 человека, должно работать хорошо в большинстве случаев.$0.267\,747\,907\,805\,267$$n, m, d.$$O(n^2\log(d))$$d$$n$

---

Объяснение решения

Функция генерации вероятности (pgf) для результатов независимых бросков кубика с стороной$n$$d$

$$d^{-n}f_n(x_1,x_2,\ldots,x_d) = d^{-n}(x_1+x_2+ \cdots + x_d)^n.$$

Коэффициент в разложении этого многочлена дает количество способов, с помощью которых могу появиться ровно раз, я е я я = 1 , 2 , ... , д $x_1^{e_1}x_2^{e_2}\cdots x_d^{e_d}$$i$$e_i$$i=1, 2, \ldots, d.$

Ограничение нашего интереса не более чем появлениями любого лица равносильно оценке по модулю идеала порожденного Чтобы выполнить эту оценку, рекурсивно используйте теорему о биномf n I x m + 1 1 , x m + 1 2 , … , x m + 1 d$m$$f_n$$\mathcal I$$x_1^{m+1}, x_2^{m+1}, \ldots, x_d^{m+1}.$

$$\eqalign{ f_n(x_1, \ldots, x_d) &= ((x_1+\cdots+x_r) + (x_{r+1}+x_{r+2} + \cdots + x_{2r}))^n \\ &= \sum_{k=0}^n \binom{n}{k} (x_1+\cdots+x_r)^k (x_{r+1}+\cdots+x_{2r})^{n-k} \\ &= \sum_{k=0}^n \binom{n}{k} f_k(x_1, \ldots, x_r) f_{n-k}(x_{r+1}, \ldots, x_{2r}) }$$

когда четно. Написание ( терминов), мы имеем$d=2r$$f_n^{(d)} = f_n(1,1,\ldots, 1)$$d$

$$f_n^{(2r)} = \sum_{k=0}^n \binom{n}{k} f_k^{(r)} f_{n-k}^{(r)}.\tag{a}$$

Если нечетно, используйте аналогичное разложение$d=2r+1$

$$\eqalign{ f_n(x_1, \ldots, x_d) &= ((x_1+\cdots+x_{2r}) + x_{2r+1})^n \\ &= \sum_{k=0}^n \binom{n}{k} f_k(x_1, \ldots, x_{2r}) f_{n-k}(x_{2r+1}), }$$

дающий

$$f_n^{(2r+1)} = \sum_{k=0}^n \binom{n}{k} f_k^{(2r)} f_{n-k}^{(1)}.\tag{b}$$

В обоих случаях мы также можем уменьшить все по модулю , что легко выполняется, начиная с$\mathcal I$

$$f_n(x_j) \cong \left\{ \matrix{x^n & n \le m \\ 0 & n \gt m} \right. \mod \mathcal{I},$$

предоставление начальных значений для рекурсии,

$$f_n^{(1)} = \left\{ \matrix{1 & n \le m \\ 0 & n \gt m} \right.$$

Что делает это эффективным, так это то, что, разбивая переменные на две равные по размеру группы переменных каждая и устанавливая все значения переменных в нам нужно только один раз оценить все для одной группы, а затем объединить результаты. Это требует вычисления до членов, каждый из которых требует вычисления для комбинации. Нам даже не нужен двумерный массив для хранения , потому что при вычислении требуются только и .$d$$r$$1,$$n+1$$O(n)$$f_n^{(r)}$$f_n^{(d)},$$f_n^{(r)}$$f_n^{(1)}$

Общее количество шагов на единицу меньше количества цифр в двоичном разложении (которое считает разбиения на равные группы в формуле ) плюс количество единиц в разложении (которое считает все разы нечетными встречается значение, требующее применения формулы ). Это все еще только шагов.$d$$(a)$$(b)$$O(\log(d))$

На Rдесятилетней рабочей станции работа была выполнена за 0,007 секунды. Код указан в конце этого поста. Он использует логарифмы вероятностей, а не сами вероятности, чтобы избежать возможных переполнений или накопления слишком больших потерь. Это позволяет удалить фактор в решении, чтобы мы могли вычислить значения, лежащие в основе вероятностей.$d^{-n}$

Обратите внимание, что эта процедура приводит к одновременному вычислению всей последовательности вероятностей , что позволяет легко изучить, как шансы изменяются с .$f_0, f_1, \ldots, f_n$$n$

---

Приложения

Распределение в обобщенной задаче дня рождения вычисляется функцией tmultinom.full. Единственная проблема заключается в поиске верхней границы для числа людей, которые должны присутствовать, прежде чем вероятность столкновения станет слишком большой. Следующий код делает это грубой силой, начиная с малого и удваивая его, пока он не станет достаточно большим. Таким образом, весь расчет занимает время где - решение. Вычисляется полное распределение вероятностей для числа людей до .$m+1$$n$$O(n^2\log(n)\log(d))$$n$$n$

#
# The birthday problem: find the number of people where the chance of
# a collision of `m+1` birthdays first exceeds `alpha`.
#
birthday <- function(m=1, d=365, alpha=0.50) {
  n <- 8
  while((p <- tmultinom.full(n, m, d))[n] > alpha) n <- n * 2
  return(p)
}

Например, минимальное количество людей, нуждающихся в толпе, чтобы было более вероятно, что по крайней мере восемь из них делят день рождения, составляет , как выяснилось из расчета . Это займет всего пару секунд. Вот график части вывода:$798$birthday(7)

---

Специальная версия этой проблемы рассматривается при распространении парадокса дня рождения более чем на 2 человека , что касается случая стороннего штампа, который бросается очень много раз.$365$

---

Код

# Compute the chance that in `n` independent rolls of a `d`-sided die, 
# no side appears more than `m` times.
#
tmultinom <- function(n, m, d, count=FALSE) tmultinom.full(n, m, d, count)[n+1]
#
# Compute the chances that in 0, 1, 2, ..., `n` independent rolls of a
# `d`-sided die, no side appears more than `m` times.
#
tmultinom.full <- function(n, m, d, count=FALSE) {
  if (n < 0) return(numeric(0))
  one <- rep(1.0, n+1); names(one) <- 0:n
  if (d <= 0 || m >= n) return(one)

  if(count) log.p <- 0 else log.p <- -log(d)
  f <- function(n, m, d) {                   # The recursive solution
    if (d==1) return(one)                    # Base case
    r <- floor(d/2)
    x <- double(f(n, m, r), m)               # Combine two equal values
    if (2*r < d) x <- combine(x, one, m)     # Treat odd `d`
    return(x)
  }
  one <- c(log.p*(0:m), rep(-Inf, n-m))      # Reduction modulo x^(m+1)
  double <- function(x, m) combine(x, x, m)
  combine <- function(x, y, m) {             # The Binomial Theorem
    z <- sapply(1:length(x), function(n) {   # Need all powers 0..n
      z <- x[1:n] + lchoose(n-1, 1:n-1) + y[n:1]
      z.max <- max(z)
      log(sum(exp(z - z.max), na.rm=TRUE)) + z.max
    })
    return(z)
  }
  x <- exp(f(n, m, d)); names(x) <- 0:n
  return(x)
}

Ответ получен с

print(tmultinom(100,20,6), digits=15)

+0,267747907805267

Метод случайной выборки

Я запускал этот код в R, повторяя 100 бросков кубика в миллион раз:

y <- репликация (1000000, все (таблица (образец (1: 6, размер = 100, замена = ИСТИНА)) <= 20))

Вывод кода внутри функции репликации - истина, если все грани появляются меньше или равны 20 раз. y - вектор с 1 миллионом значений true или false.

Общее количество истинных значений в y, деленных на 1 миллион, должны быть примерно равны вероятности, которую вы хотите. В моем случае это было 266872/1000000, что предполагает вероятность около 26,6%.

Расчет грубой силы

Этот код занимает несколько секунд на моем ноутбуке

total = 0
pb <- txtProgressBar(min = 0, max = 20^2, style = 3)
for (i in 0:20) {
  for (j in 0:20) {
    for (k in 0:20) { 
      for (l in 0:20) {
        for (m in 0:20) {
          n = 100-sum(i,j,k,l,m)
          if (n<=20) {
            total = total+dmultinom(c(i,j,k,l,m,n),100,prob=rep(1/6,6))
          }
        }
      }
    }
    setTxtProgressBar(pb, i*20+j) # update progression bar            
  }
}
total

выход: 0.2677479

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

По крайней мере, это вычисление дает упрощенно вычисленное, но действительное число для проверки других (более сложных) методов.