Estou tentando entender esse problema.
Um dado é rolado 100 vezes. Qual é a probabilidade de nenhum rosto aparecer mais de 20 vezes? Meu primeiro pensamento foi usar a distribuição binomial P (x) = 1 - 6 cmf (100, 1/6, 20), mas isso está obviamente errado, pois contamos alguns casos mais de uma vez. Minha segunda idéia é enumerar todos os rolos possíveis x1 + x2 + x3 + x4 + x5 + x6 = 100, de modo que xi <= 20 e somar os multinomiais, mas isso parece muito computacionalmente intenso. Soluções aproximadas também funcionarão para mim.

Esta é uma generalização do famoso Problema do Aniversário : dado indivíduos que têm "aniversários" aleatórios e uniformemente distribuídos entre um conjunto de possibilidades, qual é a chance de nenhum aniversário ser compartilhado por mais de indivíduos?d = 6 m = $n=100$$d=6$$m=20$

Um cálculo exato fornece a resposta (para dobrar a precisão). Vou esboçar a teoria e fornecer o código para geral O tempo assintótico do código é que o torna adequado para um número muito grande de aniversários e fornece um desempenho razoável até que esteja na casa dos milhares. Nesse ponto, a aproximação de Poisson discutida em Estendendo o paradoxo do aniversário para mais de duas pessoas deve funcionar bem na maioria dos casos.n , m , d . O ( n 2 log ( d ) ) d $0.267\,747\,907\,805\,267$$n, m, d.$$O(n^2\log(d))$$d$$n$

---

Explicação da solução

A função geradora de probabilidade (pgf) para os resultados de rolos independentes de um dado com face é$n$$d$

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

O coeficiente de na expansão desse multinomial fornece o número de maneiras pelas quais posso aparecer exatamente vezes, i e i i = 1 , 2 , ... , d $x_1^{e_1}x_2^{e_2}\cdots x_d^{e_d}$$i$$e_i$$i=1, 2, \ldots, d.$

Limitar nosso interesse a não mais que aparências por qualquer face é o mesmo que avaliar modulo o ideal gerado por Para executar esta avaliação, use o Teorema Binomial recursivamente para obterf 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}) }$$

quando é par. Escrevendo ( terms), temos$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}$$

Quando é ímpar, use uma decomposição análoga$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}), }$$

dando

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

Nos dois casos, também podemos reduzir tudo o módulo , que é facilmente realizado começando com$\mathcal I$

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

fornecendo os valores iniciais para a recursão,

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

O que torna isso eficiente é que, dividindo as variáveis em dois grupos de tamanho igual de variáveis ​​cada e definindo todos os valores das variáveis ​​como precisamos avaliar tudo apenas uma vez para um grupo e depois combinar os resultados. Isso exige a computação de até termos, cada um deles precisando de cálculo de para a combinação. Nem precisamos de uma matriz 2D para armazenar , porque ao calcular apenas e são necessários.$d$$r$$1,$$n+1$$O(n)$$f_n^{(r)}$$f_n^{(d)},$$f_n^{(r)}$$f_n^{(1)}$

O número total de etapas é um a menos que o número de dígitos na expansão binária de (que conta as divisões em grupos iguais na fórmula ) mais o número de unidades na expansão (que conta todas as vezes um número ímpar é encontrado um valor, exigindo a aplicação da fórmula ). Isso ainda é apenas etapas.$d$$(a)$$(b)$$O(\log(d))$

Em Ruma estação de trabalho de uma década, o trabalho foi realizado em 0,007 segundos. O código está listado no final desta postagem. Ele usa logaritmos das probabilidades, em vez das próprias probabilidades, para evitar possíveis estouros ou acumular muito subfluxo. Isso torna possível remover o fator na solução, para que possamos calcular as contagens subjacentes às probabilidades.$d^{-n}$

Observe que esse procedimento resulta no cálculo de toda a sequência de probabilidades uma só vez, o que facilmente nos permite estudar como as chances mudam com .$f_0, f_1, \ldots, f_n$$n$

---

Formulários

A distribuição no Problema Generalizado de Aniversário é calculada pela função tmultinom.full. O único desafio está em encontrar um limite superior para o número de pessoas que devem estar presentes antes que a chance de uma colisão se torne muito grande. O código a seguir faz isso com força bruta, começando com pequeno e dobrando-o até que seja grande o suficiente. Portanto, todo o cálculo leva tempo em que é a solução. Toda a distribuição de probabilidades para o número de pessoas até é calculada.$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)
}

Como exemplo, o número mínimo de pessoas necessárias em uma multidão para tornar mais provável que pelo menos oito delas compartilhem um aniversário é , conforme encontrado pelo cálculo . Demora apenas alguns segundos. Aqui está um gráfico de parte da saída:$798$birthday(7)

---

Uma versão especial desse problema é abordada em Estendendo o paradoxo do aniversário para mais de duas pessoas , que diz respeito ao caso de um dado de lados que é rolado muitas vezes.$365$

---

Código

# 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)
}

A resposta é obtida com

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

0.267747907805267

Método de amostragem aleatória

Eu corri esse código em R replicando 100 jogadas de dados por um milhão de vezes:

y <- replicar (1000000, todos (tabela (amostra (1: 6, tamanho = 100, substituir = VERDADEIRO)) <= 20))

A saída do código dentro da função de replicação será verdadeira se todas as faces aparecerem menores ou iguais a 20 vezes. y é um vetor com 1 milhão de valores de verdadeiro ou falso.

O total não. dos valores verdadeiros em y dividido por 1 milhão deve ser aproximadamente igual à probabilidade desejada. No meu caso, era 266872/1000000, sugerindo uma probabilidade de cerca de 26,6%

Cálculo da força bruta

Esse código leva alguns segundos no meu laptop

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

saída: 0.2677479

Ainda assim, pode ser interessante encontrar um método mais direto, caso você deseje fazer muitos desses cálculos ou usar valores mais altos, ou apenas para obter um método mais elegante.

Pelo menos esse cálculo fornece um número calculado de maneira simplista, mas válida, para verificar outros métodos (mais complicados).