J'essaie de comprendre ce problème.
Un dé est lancé 100 fois. Quelle est la probabilité qu'aucun visage n'apparaisse plus de 20 fois? Ma première pensée a été d'utiliser la distribution binomiale P (x) = 1 - 6 cmf (100, 1/6, 20) mais c'est évidemment faux puisque nous comptons certains cas plus d'une fois. Ma deuxième idée est d'énumérer tous les rouleaux possibles x1 + x2 + x3 + x4 + x5 + x6 = 100, tels que xi <= 20 et additionner les multinomiaux mais cela semble trop intensif en calcul. Des solutions approximatives fonctionneront également pour moi.

Ceci est une généralisation du fameux problème d'anniversaire : étant donné $n=100$ individus qui ont des "anniversaires" aléatoires et uniformément répartis parmi un ensemble de $d=6$ possibilités, quelle est la probabilité qu'aucun anniversaire ne soit partagé par plus de $m=20$ individus?

Un calcul exact donne la réponse (à double précision). Je vais esquisser la théorie et fournir le code général n , m , d . La synchronisation asymptotique du code est O ( n 2 log ( d ) ), ce qui le rend approprié pour un très grand nombre d'anniversaires d et fournit des performances raisonnables jusqu'à ce que n soit des milliers. À ce stade, l'approximation de Poisson discutée lors de l'extension du paradoxe d'anniversaire à plus de 2 personnesdevrait bien fonctionner dans la plupart des cas.$0.267\,747\,907\,805\,267$$n, m, d.$$O(n^2\log(d))$$d$$n$

---

Explication de la solution

La fonction de génération de probabilité (PGF) pour les résultats des bobines indépendantes de d à flancs Die$n$$d$

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

Le coefficient de dans l'expansion de ce multinomial donne le nombre de façons dont la face i peut apparaître exactement e i fois, 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.$

Limiter notre intérêt à pas plus de apparitions par n'importe quel visage revient à évaluer f n modulo l'idéal I généré par x m + 1 1 , x m + 1 2 , … , x m + 1 d . Pour effectuer cette évaluation, utilisez le théorème binomial récursivement pour obtenir$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}) }$$

lorsque est pair. En écrivant f ( d ) n = f n ( 1 , 1 , … , 1 ) ( d termes), on a$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}$$

Lorsque est impair, utilisez une décomposition analogue$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}), }$$

donnant

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

Dans les deux cas, nous pouvons également réduire tout qui est modulo I , ce qui est facile à réaliser en commençant par$\mathcal I$

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

fournir les valeurs de départ pour la récursivité,

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

Ce qui rend cela efficace, c'est qu'en divisant les variables en deux groupes de variables r de taille égale chacune et en définissant toutes les valeurs des variables à 1 , nous n'avons qu'à tout évaluer une fois pour un groupe, puis à combiner les résultats. Cela nécessite de calculer jusqu'à n + 1 termes, chacun d'eux nécessitant un calcul O ( n ) pour la combinaison. Nous n'avons même pas besoin d'un tableau 2D pour stocker le f ( r ) n , car lors du calcul de f ( d ) n , seulement $d$$r$$1,$$n+1$$O(n)$$f_n^{(r)}$$f_n^{(d)},$ etf ( 1 ) n sont requis.$f_n^{(r)}$$f_n^{(1)}$

Le nombre total d'étapes est inférieur de un au nombre de chiffres de l'expansion binaire de (qui compte les divisions en groupes égaux dans la formule ( a ) ) plus le nombre de celles de l'expansion (qui compte toutes les fois une impaire est rencontrée, nécessitant l'application de la formule ( b ) ). Ce n'est encore que des étapes O ( log ( d ) ) .$d$$(a)$$(b)$$O(\log(d))$

Sur Run poste de travail vieux de dix ans, le travail a été effectué en 0,007 secondes. Le code est répertorié à la fin de cet article. Il utilise des logarithmes des probabilités, plutôt que les probabilités elles-mêmes, pour éviter d'éventuels débordements ou accumuler trop de débordements. Cela permet de supprimer le facteur dans la solution afin que nous puissions calculer les comptes qui sous-tendent les probabilités.$d^{-n}$

Notez que cette procédure aboutit au calcul de la séquence entière des probabilités à la fois, ce qui nous permet facilement d'étudier comment les chances changent avec n .$f_0, f_1, \ldots, f_n$$n$

---

Applications

La distribution dans le problème d'anniversaire généralisé est calculée par la fonction tmultinom.full. Le seul défi consiste à trouver une limite supérieure pour le nombre de personnes qui doivent être présentes avant que les chances d'une collision ne deviennent trop grandes. Le code suivant le fait par force brute, en commençant par un petit n et en le doublant jusqu'à ce qu'il soit suffisamment grand. L'ensemble du calcul prend donc O ( n 2 log ( n ) log ( d ) ) temps où n est la solution. La distribution entière des probabilités pour le nombre de personnes jusqu'à n est calculée.$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)
}

À titre d'exemple, le nombre minimum de personnes nécessaires dans une foule pour qu'il soit plus probable qu'improbable qu'au moins huit d'entre elles partagent un anniversaire est de , comme le révèle le calcul . Cela ne prend que quelques secondes. Voici un tracé d'une partie de la sortie:$798$birthday(7)

---

Une version spéciale de ce problème est abordée dans Etendre le paradoxe de l'anniversaire à plus de 2 personnes , ce qui concerne le cas d'un dé faces qui est lancé un très grand nombre de fois.$365$

---

Code

# 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 réponse est obtenue avec

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

0,267747907805267

Méthode d'échantillonnage aléatoire

J'ai exécuté ce code dans R en répliquant 100 lancers de sorts pour un million de fois:

y <- répliquer (1000000, tout (tableau (échantillon (1: 6, taille = 100, remplacer = VRAI)) <= 20))

La sortie du code à l'intérieur de la fonction de réplication est vraie si toutes les faces apparaissent inférieures ou égales à 20 fois. y est un vecteur avec 1 million de valeurs vrai ou faux.

Le nombre total. des valeurs vraies en y divisées par 1 million devraient être approximativement égales à la probabilité que vous désirez. Dans mon cas, c'était 266872/1000000, ce qui suggère une probabilité d'environ 26,6%

Calcul de la force brute

Ce code prend quelques secondes sur mon ordinateur portable

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

sortie: 0.2677479

Mais il pourrait être intéressant de trouver une méthode plus directe au cas où vous souhaiteriez faire beaucoup de ces calculs ou utiliser des valeurs plus élevées, ou simplement pour obtenir une méthode plus élégante.

Au moins, ce calcul donne un nombre calculé de manière simpliste, mais valide, pour vérifier d'autres méthodes (plus compliquées).