擲骰子的預期次數，直到每邊出現 3 次
在每一面出現 3 次之前，您必須擲骰子的預期次數是多少？
這個問題是在新西蘭的小學提出的，並通過模擬解決。這個問題的解析解是什麼？
假設所有雙方機會均等。讓我們概括並找到直到 side 所需的預期滾動數已經出現了次，邊已經出現了次，…，和邊已經出現了次。因為雙方的身份無關緊要（他們都有相同的機會），所以可以濃縮這個目標的描述：讓我們假設側面根本不必出現，的邊只需要出現一次，…，並且雙方必須出現次。讓
指定這種情況並寫為預期的捲數。問題要求：表示所有六個面都需要被看到 3 次。 一個簡單的複發是可用的。 在下一次滾動中，出現的一側對應於其中一個: 也就是說，要么我們不需要看，要么我們需要看一次，……，或者我們需要看更多次。是我們需要看到它的次數。
- 什麼時候，我們不需要看到它，也沒有任何變化。這很有可能發生.
- 什麼時候那麼我們確實需要看到這一面。現在少了一個需要看的一面需要看到的時間和另一面次。因此，變成和變成. 讓這個操作在組件上被指定， 以便
這很有可能發生.
我們只需要計算這個骰子並使用遞歸來告訴我們還需要多少骰子。 根據期望定律和全概率定律，
（讓我們明白，無論何時，總和中的對應項為零。）
如果, 我們完成了. 否則我們可以解決，給出所需的遞歸公式
請注意
是我們希望看到的事件總數。操作對於任何一個，將該數量減少一假如，情況總是如此。因此，此遞歸在精確的深度處終止（等於在問題中）。此外（不難檢查）這個問題中每個遞歸深度的可能性數量很小（從不超過）。因此，這是一種有效的方法，至少當組合可能性不是太多並且我們記住中間結果時（因此沒有計算不止一次）。 我計算
這對我來說似乎非常小，所以我運行了一個模擬（使用
R）。擲骰數超過 300 萬次後，這款遊戲被玩到了 100,000 多次，平均時長為. 該估計的標準誤差是：該平均值與理論值的差異不大，證實了理論值的準確性。長度的分佈可能是有意義的。（顯然它必須從, 收集所有六個面所需的最小卷數，每次 3 次。）
# 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))
執行
雖然遞歸計算很簡單，它在某些計算環境中提出了一些挑戰。其中最主要的是存儲因為它們是計算出來的。這是必不可少的，否則每個值都將（冗餘地）計算很多次。但是，索引為的數組可能需要的存儲空間可能是巨大的。理想情況下，只有值在計算過程中實際遇到的應該被存儲。這需要一種關聯數組。
為了說明，這裡是工作
R代碼。註釋描述了用於存儲中間結果的簡單“AA”（關聯數組）類的創建。矢量圖被轉換為字符串，這些用於索引到E將包含所有值的列表中。這操作實現為%.%.這些預備程序啟用遞歸功能以與數學符號平行的方式相當簡單地定義。特別是，線
x <- (d + sum(sapply(1:n, function(i) j[i+1]*e.(j %.% i))))/(d - j[1])與公式直接可比多於。請注意，所有索引都增加了因為
R開始索引它的數組在而不是.時間表明它需要計算秒數
e(c(0,0,0,6))；它的價值是32.6771634160506
累積的浮點舍入誤差破壞了最後兩位數字（應該是
68而不是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))最後，這是產生準確答案的原始Mathematica實現。記憶是通過慣用的
e[i_] := e[i] = ...表達方式完成的，幾乎消除了所有的R預備。但是，在內部，這兩個程序以相同的方式做同樣的事情。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}]
