Estoy tratando de envolver mi cabeza alrededor de este problema.
Un dado es lanzado 100 veces. ¿Cuál es la probabilidad de que no se aparezca más de 20 veces? Mi primer pensamiento fue el uso de la distribución Binomial P(x) = 1 - 6 cmf(100, 1/6, 20), pero esto es evidentemente falso, ya que el recuento de algunos casos más de una vez. Mi segunda idea es la de enumerar todos los posibles rollos de x1 +x2+x3+x4+x5+x6 = 100, tal que xi <= 20 y la suma de la multinomials pero esto parece demasiado intensas. Aproximado soluciones también han de trabajar para mí.
Respuestas
¿Demasiados anuncios?
jldugger
Puntos
7490
Esta es una generalización de la famosa Cumpleaños Problema: dado $n=100$ de los individuos que han aleatorios uniformemente distribuidos "cumpleaños" entre un conjunto de $d=6$ posibilidades, ¿cuál es la probabilidad de que ninguna de cumpleaños es compartida por más de $m=20$ de los individuos?
Un cálculo exacto de los rendimientos de la respuesta $0.267\,747\,907\,805\,267$ (doble precisión). Voy a esbozar la teoría y proporcionar el código general $n, m, d.$ El código del asintótica tiempo es $O(n^2\log(d))$ que lo hace adecuado para grandes cantidades de cumpleaños $d$ y proporciona razonable de rendimiento de hasta el $n$ está en los miles. En ese momento, la aproximación de Poisson discutido en la Ampliación de la paradoja de cumpleaños para más de 2 personas debería funcionar bien en la mayoría de los casos.
Explicación de la solución
La probabilidad de generación de función (pgf) para los resultados de la $n$ independiente rollos de una $d$colindado mueren es
$$d^{-n}f_n(x_1,x_2,\ldots,x_d) = d^{-n}(x_1+x_2+ \cdots + x_d)^n.$$
The coefficient of $x_1^{e_1}x_2^{e_2}\cdots x_d^{e_d}$ in the expansion of this multinomial gives the number of ways in which face $i$ can appear exactly $e_i$ times, $i=1, 2, \ldots, d.$
Limiting our interest to no more than $m$ appearances by any face is tantamount to evaluating $f_n$ modulo the ideal $\mathcal I$ generated by $x_1^{m+1}, x_2^{m+1}, \ldots, x_d^{m+1}.$ To perform this evaluation, use the Binomial Theorem recursively to obtain
$$\eqalign{ f_n(x_1, \ldots, x_d) y= ((x_1+\cdots+x_r) + (x_{i+1}+x_{i+2} + \cdots + x_{2r}))^n \\ &= \sum_{k=0}^n \binom{n}{k} (x_1+\cdots+x_r)^k (x_{i+1}+\cdots+x_{2r})^{n-k} \\ &= \sum_{k=0}^n \binom{n}{k} f_k(x_1, \ldots, x_r) f_{n-k}(x_{i+1}, \ldots, x_{2r}) }$$
when $d=2r$ is even. Writing $f_n^{(d)} = f_n(1,1,\ldots, 1)$ ($d$ terms), we have
$$f_n^{(2r)} = \sum_{k=0}^n \binom{n}{k} f_k^{(r)} f_{n-k}^{(r)}.\tag{a}$$
When $d=2r+1$ is odd, use an analogous decomposition
$$\eqalign{ f_n(x_1, \ldots, x_d) y= ((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}), }$$
giving
$$f_n^{(2r+1)} = \sum_{k=0}^n \binom{n}{k} f_k^{(2r)} f_{n-k}^{(1)}.\tag{b}$$
In both cases, we may also reduce everything modulo $\mathcal I$, which is easily carried out beginning with
$$f_n(x_j) \cong \left\{ \matrix{x^n & n \le m \\ 0 & n \gt m} \right. \mod \mathcal{I},$$
providing the starting values for the recursion,
$$f_n^{(1)} = \left\{ \matrix{1 & n \le m \\ 0 & n \gt m} \right.$$
What makes this efficient is that by splitting the $d$ variables into two equal-sized groups of $r$ variables each and setting all variable values to $1,$ we only have to evaluate everything once for one group and then combine the results. This requires computing up to $n+1$ terms, each of them needing $O(n)$ calculation for the combination. We don't even need a 2D array to store the $f_n^{(r)}$, because when computing $f_n^{(d)},$ only $f_n^{(r)}$ and $f_n^{(1)}$ are required.
The total number of steps is one less than the number of digits in the binary expansion of $d$ (which counts the splits into equal groups in formula $(a)$) plus the number of ones in the expansion (which counts all of the times an odd value is encountered, requiring the application of formula $(b)$). That's still just $O(\log(d))$ steps.
In R on a decade-old workstation the work was done in 0.007 seconds. The code is listed at the end of this post. It uses logarithms of the probabilities, rather than the probabilities themselves, to avoid possible overflows or accumulating too much underflow. This makes it possible to remove the $d^{-n}$ factor in the solution so we may compute the counts that underlie the probabilities.
Note that this procedure results in computing the whole sequence of probabilities $f_0, f_1, \ldots, f_n$ at once, which easily enables us to study how the chances change with $n$.
Applications
The distribution in the generalized Birthday Problem is computed by the function tmultinom.full. The only challenge lies in finding an upper bound for the number of people who must be present before the chance of an $m+1$-collision becomes too great. The following code does this by brute force, starting with small $n$ and doubling it until it's large enough. The whole calculation therefore takes $O(n^2\log(n)\log(d))$ time where $n$ is the solution. The entire distribution of probabilities for numbers of people up through $n$ is computed.
#
# 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)
}
As an example, the minimum number of people needed in a crowd to make it more likely than not that at least eight of them share a birthday is $798$, as found by the calculation birthday(7). It takes just a couple of seconds. Here is a plot of part of the output:
A special version of this problem is addressed at Extending the birthday paradox to more than 2 people, which concerns the case of a $365$-cara de morir, que se rodó un número muy grande de veces.
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
user164061
Puntos
281
La fuerza bruta de cálculo
Este código se toma un par de segundos en mi 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
salida: 0.2677479
Pero todavía podría ser interesante para encontrar un método más directo en caso de que usted desee hacer un montón de estos cálculos o utilice valores más altos, o simplemente por el bien de conseguir una más elegante método.
Al menos este cálculo da un resumen calculado, pero válido, el número de consultar otras (más complicado) métodos.
Método de muestreo aleatorio
Me encontré este código en R replicando 100 mueren a tiros por un millón de veces:
y <- replicate(1000000, todos(de la tabla(de la muestra(1:6, size = 100, replace = TRUE)) <=20))
La salida del código dentro de la función replicate es verdadera si todas las caras que aparecen inferior o igual a 20 veces. y es un vector con 1 millón de valores de verdadero o falso.
El total no. de los verdaderos valores en y dividido por 1 millón debe ser aproximadamente igual a la probabilidad de que el deseo. En mi caso fue 266872/1000000, lo que sugiere una probabilidad de alrededor de 26.6%
