Próbuję obejść ten problem.
Kostka jest rzucana 100 razy. Jakie jest prawdopodobieństwo, że żadna twarz nie pojawi się więcej niż 20 razy? Moją pierwszą myślą było użycie rozkładu dwumianowego P (x) = 1 - 6 cmf (100, 1/6, 20), ale jest to oczywiście błędne, ponieważ liczymy niektóre przypadki więcej niż jeden raz. Moim drugim pomysłem jest wyliczenie wszystkich możliwych rolek x1 + x2 + x3 + x4 + x5 + x6 = 100, tak aby xi <= 20 i zsumować wielomiany, ale wydaje się to zbyt intensywne obliczeniowo. Przybliżone rozwiązania również będą dla mnie działać.

Jest to uogólnienie znanego problemu urodzinowego : biorąc pod uwagę $n=100$ osób, które mają losowe, równomiernie rozmieszczone „urodziny” wśród zestawu $d=6$ możliwości, jaka jest szansa, że ​​żadne urodziny nie są dzielone przez więcej niż $m=20$ osób?

Dokładne obliczenia dają odpowiedź (z podwójną precyzją). Naszkicuję teorię i podam kod dla ogólnych n , m , d . Czas asymptotyczny kodu wynosi O ( n 2 log ( d ) ), co czyni go odpowiednim dla bardzo dużej liczby urodzin d i zapewnia rozsądną wydajność, dopóki n nie będzie tysiąca. W tym momencie aproksymacja Poissona omawiana podczasrozszerzania paradoksu urodzinowego na więcej niż 2 osobypowinna w większości przypadków działać dobrze.$0.267\,747\,907\,805\,267$$n, m, d.$$O(n^2\log(d))$$d$$n$

---

Objaśnienie rozwiązania

Funkcja generowania prawdopodobieństwa (PGF) w odniesieniu do wyników niezależne rolki o d -sided dyszowej mieści$n$$d$

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

Współczynnik w rozwinięciu tego wielomianu daje liczbę sposobów, w jakie twarz i może pojawić się dokładnie e i razy, 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.$

Ograniczenie naszego zainteresowania nie więcej niż występy dowolnej powierzchni jest równoznaczne oceniającej f n modulo idealne I generowany przez x m + 1 1 , x m + 1, 2 , ... , x m + 1 d . Aby wykonać tę ocenę, użyj rekurencyjnie twierdzenia dwumianowego w celu uzyskania$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}) }$$

gdy jest parzyste. Pisząc ( terminy), mamyf ( d ) n = f n ( 1 , 1 , … , 1 ) $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}$$

Gdy jest nieparzyste, zastosuj analogiczny rozkład$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}), }$$

dający

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

W obu przypadkach możemy również zredukować wszystko modulo , co można łatwo przeprowadzić od początku$\mathcal I$

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

podając wartości początkowe rekurencji,

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

To, co sprawia, że ​​jest to wydajne, polega na tym, że dzieląc zmiennych na dwie równe grupy zmiennych każda i ustawiając wszystkie wartości zmiennych na musimy ocenić wszystko tylko raz dla jednej grupy, a następnie połączyć wyniki. Wymaga to obliczenia do terminów, z których każdy wymaga obliczenia dla kombinacji. Nie potrzebujemy nawet tablicy 2D do przechowywania , ponieważ podczas obliczania wymagane są tylko i .r 1 , n + 1 O ( n ) f ( r ) n f ( d ) n , f ( r ) n f ( 1 )$d$$r$$1,$$n+1$$O(n)$$f_n^{(r)}$$f_n^{(d)},$$f_n^{(r)}$$f_n^{(1)}$

Całkowita liczba kroków jest o jeden mniejsza niż liczba cyfr w binarnym rozwinięciu (który liczy podziały na równe grupy we wzorze ) plus liczba jednych w rozwinięciu (które liczą wszystkie razy nieparzyste napotkano wartość wymagającą zastosowania wzoru ). To wciąż tylko kroki .( a ) ( b ) O ( log ( d ) $d$$(a)$$(b)$$O(\log(d))$

Na Rdziesięcioletniej stacji roboczej praca została wykonana w 0,007 sekundy. Kod znajduje się na końcu tego postu. Wykorzystuje logarytmy prawdopodobieństwa, a nie same prawdopodobieństwa, aby uniknąć ewentualnych przelewów lub nagromadzenia zbyt dużego niedomiaru. Umożliwia to usunięcie czynnika w rozwiązaniu, dzięki czemu możemy obliczyć liczby leżące u podstaw prawdopodobieństw.$d^{-n}$

Zauważ, że ta procedura skutkuje obliczeniem całej sekwencji prawdopodobieństwa jednocześnie, co pozwala nam łatwo zbadać, jak szanse zmieniają się z .$f_0, f_1, \ldots, f_n$$n$

---

Aplikacje

Rozkład w uogólnionym problemie urodzinowym jest obliczany przez funkcję tmultinom.full. Jedyne wyzwanie polega na znalezieniu górnej granicy dla liczby osób, które muszą być obecne, zanim szansa kolizji stanie się zbyt duża. Poniższy kod robi to brutalną siłą, zaczynając od małego i podwajając go, aż będzie wystarczająco duży. Całe obliczenie zajmuje zatem czas gdzie jest rozwiązaniem. Obliczany jest cały rozkład prawdopodobieństwa dla liczby osób w górę przez .n O ( n 2 log ( n ) log ( d ) ) n $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)
}

Przykładowo, jak wynika z obliczeń , minimalna liczba osób potrzebnych w tłumie, aby zwiększyć prawdopodobieństwo, że co najmniej osiem z nich ma wspólne urodziny, wynosi . To zajmuje tylko kilka sekund. Oto wykres części wyniku:$798$birthday(7)

---

Specjalna wersja tego problemu dotyczy rozszerzenia paradoksu urodzinowego na więcej niż 2 osoby , który dotyczy przypadku stronnej kostki, która jest rzucana bardzo wiele razy.$365$

---

Kod

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

Odpowiedź uzyskuje się za pomocą

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

0,267747907805267

Metoda losowego próbkowania

Uruchomiłem ten kod w R replikując 100 rzutów kostką milion razy:

y <- replikacja (1000000, wszystkie (tabela (próbka (1: 6, rozmiar = 100, zamień = PRAWDA)) <= 20))

Dane wyjściowe kodu w funkcji replikacji są prawdziwe, jeśli wszystkie twarze wydają się mniejsze lub równe 20 razy. y jest wektorem o 1 milionie wartości true lub false.

Łączna liczba prawdziwych wartości y podzielonych przez 1 milion powinno być w przybliżeniu równe pożądanemu prawdopodobieństwu. W moim przypadku było to 266872/1000000, co sugeruje prawdopodobieństwo około 26,6%

Obliczanie siły brutalnej

Ten kod zajmuje kilka sekund na moim laptopie

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

wyjście: 0,2677479

Ale nadal może być interesujące znalezienie bardziej bezpośredniej metody na wypadek, gdybyś chciał wykonać wiele z tych obliczeń lub użyć wyższych wartości, lub po prostu w celu uzyskania bardziej eleganckiej metody.

Przynajmniej to obliczenie daje uproszczoną, ale prawidłową liczbę do sprawdzenia innych (bardziej skomplikowanych) metod.