Estoy tratando de entender este problema.
Un dado se tira 100 veces. ¿Cuál es la probabilidad de que ninguna cara aparezca más de 20 veces? Mi primer pensamiento fue usar la distribución binomial P (x) = 1 - 6 cmf (100, 1/6, 20) pero esto obviamente es incorrecto ya que contamos algunos casos más de una vez. Mi segunda idea es enumerar todos los rollos posibles x1 + x2 + x3 + x4 + x5 + x6 = 100, de modo que xi <= 20 y sume los multinomiales, pero esto parece demasiado computacionalmente intensivo. Las soluciones aproximadas también funcionarán para mí.

Esta es una generalización del famoso problema de cumpleaños : dado individuos que tienen "cumpleaños" aleatorios y distribuidos de manera uniforme entre un conjunto de posibilidades, ¿cuál es la posibilidad de que más de individuos no compartan cumpleaños ?$n=100$$d=6$$m=20$

Un cálculo exacto arroja la respuesta (con doble precisión). Dibujaré la teoría y proporcionaré el código para El tiempo asintótico del código es lo que lo hace adecuado para un gran número de cumpleaños y proporciona un rendimiento razonable hasta que esté en los miles. En ese punto, la aproximación de Poisson discutida en Extender la paradoja del cumpleaños a más de 2 personas debería funcionar bien en la mayoría de los casos.$0.267\,747\,907\,805\,267$$n, m, d.$$O(n^2\log(d))$$d$$n$

---

Explicación de la solución.

La función de generación de probabilidad (pgf) para los resultados de tiradas independientes de un dado de lado es$n$$d$

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

El coeficiente de en la expansión de este multinomial da el número de formas en las que se enfrentan puedo aparecer exactamente veces, 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 nuestro interés a no más de apariencias por cualquier cara equivale a evaluar módulo la ideal que mediante Para realizar esta evaluación, utilice el teorema binomial de forma recursiva para obtenerf 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}) }$$

cuando es par. Escribiendo ( términos), tenemos$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}$$

Cuando es impar, use una descomposición 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}$$

En ambos casos, también podemos reducir todo el módulo , que se lleva a cabo fácilmente comenzando con$\mathcal I$

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

proporcionando los valores iniciales para la recursividad,

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

Lo que hace que esto sea eficiente es que al dividir las variables en dos grupos de variables de igual tamaño y establecer todos los valores de las variables en solo tenemos que evaluar todo una vez para un grupo y luego combinar los resultados. Esto requiere calcular hasta términos, cada uno de ellos necesita un cálculo de para la combinación. Ni siquiera necesitamos una matriz 2D para almacenar , porque al calcular solo se requieren y .$d$$r$$1,$$n+1$$O(n)$$f_n^{(r)}$$f_n^{(d)},$$f_n^{(r)}$$f_n^{(1)}$

El número total de pasos es uno menos que el número de dígitos en la expansión binaria de (que cuenta las divisiones en grupos iguales en la fórmula ) más el número de unos en la expansión (que cuenta todas las veces un número impar se encuentra un valor que requiere la aplicación de la fórmula ). Eso sigue siendo solo pasos .$d$$(a)$$(b)$$O(\log(d))$

En Runa estación de trabajo de una década, el trabajo se realizó en 0.007 segundos. El código aparece al final de esta publicación. Utiliza logaritmos de las probabilidades, en lugar de las probabilidades en sí, para evitar posibles desbordamientos o acumular demasiados desbordamientos. Esto hace posible eliminar el factor en la solución para que podamos calcular los recuentos que subyacen a las probabilidades.$d^{-n}$

Tenga en cuenta que este procedimiento da como resultado el cálculo de toda la secuencia de probabilidades a la vez, lo que nos permite estudiar fácilmente cómo cambian las posibilidades con .$f_0, f_1, \ldots, f_n$$n$

---

Aplicaciones

La distribución en el problema de cumpleaños generalizado es calculada por la función tmultinom.full. El único desafío radica en encontrar un límite superior para el número de personas que deben estar presentes antes de que la posibilidad de una colisión sea ​​demasiado grande. El siguiente código hace esto por fuerza bruta, comenzando con una pequeña y duplicándola hasta que sea lo suficientemente grande. Por lo tanto, todo el cálculo lleva tiempo donde es la solución. Se calcula la distribución completa de probabilidades para el número de personas hasta .$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 ejemplo, el número mínimo de personas que se necesitan en una multitud para que sea más probable que al menos ocho de ellos compartan un cumpleaños es , según el cálculo . Solo toma un par de segundos. Aquí hay una gráfica de parte de la salida:$798$birthday(7)

---

Una versión especial de este problema se aborda en Extender la paradoja del cumpleaños a más de 2 personas , que se refiere al caso de un dado de lados que se lanza una gran cantidad de veces.$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)
}

La respuesta se obtiene con

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

0.267747907805267

Método de muestreo aleatorio

Ejecuté este código en R replicando 100 lanzamientos de dados por un millón de veces:

y <- replicar (1000000, all (tabla (muestra (1: 6, tamaño = 100, reemplazar = VERDADERO)) <= 20))

La salida del código dentro de la función de replicación es verdadera si todas las caras aparecen menos o igual que 20 veces. y es un vector con 1 millón de valores de verdadero o falso.

El total no. de los valores verdaderos en y dividido por 1 millón debería ser aproximadamente igual a la probabilidad que desea. En mi caso fue 266872/1000000, lo que sugiere una probabilidad de alrededor del 26,6%

Cálculo de fuerza bruta

Este código tarda unos segundos en mi computadora portátil

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

salida: 0.2677479

Pero aún así puede ser interesante encontrar un método más directo en caso de que desee hacer muchos de estos cálculos o usar valores más altos, o simplemente por el simple hecho de obtener un método más elegante.

Al menos este cálculo proporciona un número calculado de manera simplista, pero válido, para verificar otros métodos (más complicados).