Oczekiwano, ile razy rzucisz kością, aż każda strona pojawi się 3 razy

Jakiej oczekiwanej liczby razy musisz rzucić kostką, aż każda ze stron pojawi się 3 razy?

To pytanie zostało zadane w szkole podstawowej w Nowej Zelandii i zostało rozwiązane za pomocą symulacji. Jakie jest analityczne rozwiązanie tego problemu?

Załóżmy, że wszystkie stron mają równe szanse. oczekiwaną liczbę potrzebnych rzutów, aż strona się razy, strona się razy, ..., a strona się razy. Ponieważ tożsamość stron nie ma znaczenia (wszystkie mają równe szanse), opis tego celu można : załóżmy, że strony wcale nie muszą się pojawiać, strony muszą pojawić się tylko raz , ... i boków muszą pojawić się razy. Pozwolić$d=6$$1$$n_1$$2$$n_2$$d$$n_d$$i_0$$i_1$$i_n$$n=\max(n_1,n_2,\ldots,n_d)$e ( i ) e ( 0 , 0 , 0 , 6

$$\mathbf{i}=(i_0,i_1,\ldots,i_n)$$ tę sytuację i napisz dla oczekiwanej liczby rzutów. Pytanie wymaga : oznacza, że ​​wszystkie sześć stron trzeba zobaczyć trzy razy. $$e(\mathbf{i})$$ $e(0,0,0,6)$$i_3 = 6$

Dostępna jest łatwa rekurencja. Na następnej stronie rolki, które się pojawi, odpowiada jeden z : To znaczy, że albo nie trzeba go zobaczyć, czy potrzebujemy go zobaczyć raz, ..., lub co potrzebne, aby go zobaczyć razy więcej . to ile razy potrzebowaliśmy to zobaczyć. n $i_j$$n$$j$

• Gdy , nie musieliśmy tego widzieć i nic się nie zmienia. Dzieje się tak z prawdopodobieństwem .i 0 /$j=0$$i_0/d$

• Kiedy musieliśmy zobaczyć tę stronę. Teraz jest o jedną stronę mniejszą niż raz i jeszcze jedną stronę, którą trzeba zobaczyć razy. Zatem staje się a staje się . Niech ta operacja na komponentach zostanie oznaczona , tak abyj j - 1 i j i j - 1 i j - 1 i j + 1 i i ⋅ $j \gt 0$$j$$j-1$$i_j$$i_j-1$$i_{j-1}$$i_j+1$$\mathbf{i}$$\mathbf{i}\cdot j$

  $$\mathbf{i}\cdot j = (\color{gray}{i_0, \ldots, i_{j-2}}, i_{j-1}+1, i_j-1, \color{gray}{i_{j+1},\ldots, i_n}).$$

  Dzieje się tak z prawdopodobieństwem .$i_j/d$

Musimy tylko policzyć ten rzut kości i użyć rekurencji, aby powiedzieć nam, ile oczekiwanych jest rzutów. Według praw oczekiwania i całkowitego prawdopodobieństwa,

$$e(\mathbf{i}) = 1 + \frac{i_0}{d}e(\mathbf{i}) + \sum_{j=1}^n \frac{i_j}{d}e(\mathbf{i}\cdot j)$$

(Rozumiemy, że ilekroć , odpowiedni warunek w sumie wynosi zero).$i_j=0$

Jeśli , to skończymy, a . W przeciwnym razie możemy rozwiązać dla , podając żądaną formułę rekurencyjnąe ( i ) = 0 e ( i$i_0=d$$e(\mathbf{i}) =0$$e(\mathbf{i})$

$$e(\mathbf{i}) = \frac{d + i_1 e(\mathbf{i}\cdot 1) + \cdots + i_n e(\mathbf{i}\cdot n)}{d - i_0}.\tag{1}$$

Zauważ, że to łączna liczba zdarzeń, które chcemy zobaczyć. Operacja zmniejsza tę liczbę o jeden dla dowolnego pod warunkiem, że , co zawsze ma miejsce. Dlatego ta rekurencja kończy się na głębokości dokładnie(równa w pytaniu). Ponadto (co nie jest trudne do sprawdzenia) liczba możliwości na każdej głębokości rekurencji w tym pytaniu jest niewielka (nigdy nie przekracza ). W związku z tym jest to skuteczna metoda, przynajmniej wtedy, gdy możliwości kombinatoryczne nie są zbyt liczne, a my zapamiętujemy wyniki pośrednie (tak aby żadna wartość⋅ j j > 0 i j > 0 | i | 3 ( 6 ) = 18 8

$$|\mathbf{i}| = 0(i_0) + 1(i_1) + \cdots + n(i_n)$$ $\cdot j$$j\gt 0$$i_j \gt 0$$|\mathbf{i}|$$3(6) = 18$$8$$e$ jest obliczany więcej niż jeden raz).

Obliczam, że

$$e(0,0,0,6) = \frac{2\,286\,878\,604\,508\,883}{69\,984\,000\,000\,000}\approx 32.677.$$

Wydawało mi się to okropnie małe, więc przeprowadziłem symulację (używając R). Po ponad trzech milionach rzutów kostką gra została ukończona ponad 100 000 razy, przy średniej długości . Standardowy błąd tego oszacowania wynosi : różnica między tą średnią a wartością teoretyczną jest nieznaczna, co potwierdza dokładność wartości teoretycznej.$32.669$$0.027$

Interesujący może być rozkład długości. (Oczywiście musi zaczynać się od , minimalna liczba rzutów potrzebna do zebrania wszystkich sześciu stron trzy razy każda).$18$

# Specify the problem
d <- 6   # Number of faces
k <- 3   # Number of times to see each
N <- 3.26772e6 # Number of rolls

# Simulate many rolls
set.seed(17)
x <- sample(1:d, N, replace=TRUE)

# Use these rolls to play the game repeatedly.
totals <- sapply(1:d, function(i) cumsum(x==i))
n <- 0
base <- rep(0, d)
i.last <- 0
n.list <- list()
for (i in 1:N) {
  if (min(totals[i, ] - base) >= k) {
    base <- totals[i, ]
    n <- n+1
    n.list[[n]] <- i - i.last
    i.last <- i
  }
}

# Summarize the results
sim <- unlist(n.list)
mean(sim)
sd(sim) / sqrt(length(sim))
length(sim)
hist(sim, main="Simulation results", xlab="Number of rolls", freq=FALSE, breaks=0:max(sim))

---

Realizacja

Chociaż rekurencyjne obliczanie jest proste, stanowi pewne wyzwanie w niektórych środowiskach obliczeniowych. Najważniejszym z nich jest przechowywanie wartości podczas ich obliczania. Jest to niezbędne, ponieważ w przeciwnym razie każda wartość zostanie (nadmiarowo) obliczona bardzo wiele razy. Jednak pamięć potencjalnie potrzebna dla tablicy indeksowanej przez może być ogromna. Idealnie powinny być przechowywane tylko wartości , które faktycznie występują podczas obliczeń. Wymaga to pewnego rodzaju tablicy asocjacyjnej.$e$$e(\mathbf{i})$$\mathbf{i}$$\mathbf{i}$

Aby to zilustrować, oto działający Rkod. Komentarze opisują utworzenie prostej klasy „AA” (tablica asocjacyjna) do przechowywania wyników pośrednich. Wektory są konwertowane na ciągi, które służą do indeksowania do listy, która będzie zawierać wszystkie wartości. Operacja jest zaimplementowana jako .$\mathbf{i}$E$\mathbf{i}\cdot j$%.%

Te czynności wstępne umożliwiają raczej zdefiniowanie funkcji rekurencyjnej w sposób podobny do zapisu matematycznego. W szczególności linia$e$

x <- (d + sum(sapply(1:n, function(i) j[i+1]*e.(j %.% i))))/(d - j[1])

jest bezpośrednio porównywalny z powyższym wzorem . Zauważ, że wszystkie indeksy zostały zwiększone o ponieważ zaczyna indeksować swoje tablice od zamiast .$(1)$$1$R$1$$0$

Czas pokazuje, że obliczenie zajmuje sekundy ; jego wartość to$0.01$e(c(0,0,0,6))

32,6771634160506

Skumulowany błąd zaokrąglenia zmiennoprzecinkowego zniszczył dwie ostatnie cyfry (które powinny być 68raczej niż 06).

e <- function(i) {
  #
  # Create a data structure to "memoize" the values.
  #
  `[[<-.AA` <- function(x, i, value) {
    class(x) <- NULL
    x[[paste(i, collapse=",")]] <- value
    class(x) <- "AA"
    x
  }
  `[[.AA` <- function(x, i) {
    class(x) <- NULL
    x[[paste(i, collapse=",")]]
  }
  E <- list()
  class(E) <- "AA"
  #
  # Define the "." operation.
  #
  `%.%` <- function(i, j) {
    i[j+1] <- i[j+1]-1
    i[j] <- i[j] + 1
    return(i)
  }
  #
  # Define a recursive version of this function.
  #
  e. <- function(j) {
    #
    # Detect initial conditions and return initial values.
    #
    if (min(j) < 0 || sum(j[-1])==0) return(0)
    #
    # Look up the value (if it has already been computed).
    #
    x <- E[[j]]
    if (!is.null(x)) return(x)
    #
    # Compute the value (for the first and only time).
    #
    d <- sum(j)
    n <- length(j) - 1
    x <- (d + sum(sapply(1:n, function(i) j[i+1]*e.(j %.% i))))/(d - j[1])
    #
    # Store the value for later re-use.
    #
    E[[j]] <<- x
    return(x)
  }
  #
  # Do the calculation.
  #
  e.(i)
}
e(c(0,0,0,6))

Wreszcie, oto oryginalna implementacja Mathematica, która dała dokładną odpowiedź. Zapamiętywanie odbywa się za pomocą e[i_] := e[i] = ...wyrażenia idiomatycznego , eliminując prawie wszystkie czynności Rwstępne. Jednak wewnętrznie oba programy robią te same rzeczy w ten sam sposób.

shift[j_, x_List] /; Length[x] >= j >= 2 := Module[{i = x},
   i[[j - 1]] = i[[j - 1]] + 1;
   i[[j]] = i[[j]] - 1;
   i];
e[i_] := e[i] = With[{i0 = First@i, d = Plus @@ i},
    (d + Sum[If[i[[k]] > 0, i[[k]]  e[shift[k, i]], 0], {k, 2, Length[i]}])/(d - i0)];
e[{x_, y__}] /; Plus[y] == 0  := e[{x, y}] = 0

e[{0, 0, 0, 6}]

$\frac{2286878604508883}{69984000000000}$

Oryginalna wersja tego pytania rozpoczęła życie od zadania:

ile rzutów potrzeba, aż każda strona pojawi się 3 razy

Oczywiście jest to pytanie, na które nie ma odpowiedzi, jak to skomentował @JuhoKokkala: odpowiedź jest zmienną losową z rozkładem, który należy znaleźć. Pytanie zostało następnie zmodyfikowane, aby zapytać: „Jaka jest oczekiwana liczba rzutów”. Poniższa odpowiedź ma na celu odpowiedzieć na postawione pierwotne pytanie: jak znaleźć rozkład liczby rolek bez korzystania z symulacji i po prostu przy użyciu prostych koncepcyjnie technik, które każdy student z Nowej Zelandii za pomocą komputera mógłby wdrożyć prawo Dlaczego nie? Problem sprowadza się do 1-liniowej.$\rightarrow$

Rozkład liczby wymaganych rolek ... tak, że każda strona pojawia się 3 razy

Rzucamy kostką razy. Niech oznacza, ile razy pojawia się strona matrycy, gdzie . Zatem wspólnym pmf dla jest tj . :$n$$X_i$$i$$i \in \{1, \dots, 6\}$$(X_1, X_2,\dots, X_6)$$\text{Multinomial}(n,\frac16)$

$$P\left(X_1=x_1,\ldots ,X_6=x_6\right) \; = \; \frac{n! }{ x_1! \cdots x_6!} \; \frac{1}{6^n} \quad \text{ subject to: } \quad \sum _{i=1}^6 x_i=n$$

Niech:Zatem cdf z to:$\quad N = \min\big\{n: \; {X_i \geq 3 \; \forall_i } \big\}. \;$$N$$\quad P(N \leq n) \; = \; P\big(X_{\forall_i} \geq 3 \; \big| \; n\big)$

tzn. Aby znaleźć cdf , po prostu oblicz dla każdej wartości :$P(N \leq n)$$n = \{18, 19, 20,\dots\}$

$$P(X_1 \geq3, \dots , X_6 \geq 3) \quad \text{ where } \quad (X_1, \dots, X_6) \sim \text{Multinomial}(n,\frac16)$$

Oto na przykład kod Mathematica , który to robi, gdy wzrasta z 18 do 60. Jest to w zasadzie jednowierszowy:$n$

 cdf = ParallelTable[ 
   Probability[x1 >= 3 && x2 >= 3 && x3 >= 3 && x4 >= 3 && x5 >= 3 &&  x6 >= 3, 
       {x1, x2, x3, x4, x5, x6} \[Distributed] MultinomialDistribution[n, Table[1/6, 6]]],
    {n, 18, 60}]

... co daje dokładny plik cdf w miarę wzrostu :$n$

$$\begin{array}{cc} 18 & \frac{14889875}{11019960576} \\ 19 & \frac{282907625}{44079842304} \\ 20 & \frac{3111983875}{176319369216} \\ 21 & \frac{116840849125}{3173748645888} \\ 22 & \frac{3283142988125}{50779978334208} \\ 23 & \frac{61483465418375}{609359740010496} \\ \vdots & \vdots\\ \\ \end{array}$$

Oto wykres cdf , w funkcji :$P(N\leq n)$$n$

Aby uzyskać pmf , po prostu najpierw rozróżnij plik cdf:$P(N=n)$

Oczywiście rozkład nie ma górnej granicy, ale możemy tutaj łatwo rozwiązać tyle wartości, ile jest praktycznie wymagane. Podejście to jest ogólne i powinno działać równie dobrze dla każdej wymaganej kombinacji wymaganych stron.