模具 100 卷 沒有臉出現超過 20 次
我正試圖解決這個問題。
一個骰子被擲了 100 次。沒有面孔出現超過 20 次的概率是多少？我的第一個想法是使用二項分佈 P(x) = 1 - 6 cmf(100, 1/6, 20) 但這顯然是錯誤的，因為我們不止一次計算某些情況。我的第二個想法是枚舉所有可能的滾動 x1 +x2+x3+x4+x5+x6 = 100，使得 xi <= 20 並對多項式求和，但這似乎計算量太大。近似的解決方案也對我有用。
這是著名的生日問題的概括：給定在一組人群中具有隨機、均勻分佈的“生日”的個人可能性，不超過生日的概率是多少個人？
精確計算得出答案（雙精度）。我將概述理論並提供通用代碼 代碼的漸近時間是這使得它適合大量的生日並提供合理的性能，直到數以千計。那時，將生日悖論擴展到超過 2 人中討論的泊松近似在大多數情況下應該可以很好地工作。
解決方案說明
結果的概率生成函數 (pgf)一個獨立的捲雙面模具是
係數在這個多項式的展開中給出了面對的方式的數量可以準確出現次，
將我們的興趣限制在不超過任何一張臉都等於評價以理想為模由產生 要執行此評估，請遞歸使用二項式定理以獲得
什麼時候甚至。寫作(條款），我們有
什麼時候是奇數，使用類似的分解
給予
在這兩種情況下，我們也可以減少一切模, 這很容易從
提供遞歸的起始值，
使這種方法有效的原因在於，通過拆分變量分成兩組大小相等的每個變量並將所有變量值設置為我們只需要為一組評估一次，然後合併結果。這需要計算高達條款，他們每個人都需要組合計算。我們甚至不需要二維數組來存儲, 因為當計算只要和是必須的。
總步數比二進制展開的位數少一（它在公式中將拆分計數為相等的組) 加上展開式中的個數（計算遇到奇數的所有次數，需要應用公式）。那還只是腳步。
在
R一個有十年曆史的工作站上，這項工作在 0.007 秒內完成。代碼列在本文末尾。它使用概率的對數，而不是概率本身，以避免可能的上溢或累積過多的下溢。這使得可以刪除解決方案中的因素，因此我們可以計算構成概率的計數。請注意，此過程會導致計算整個概率序列一次，這使我們很容易研究機會如何隨著.
應用
廣義生日問題中的分佈由函數 計算
tmultinom.full。唯一的挑戰在於找到一個上限，即在有機會出現之前必須在場的人數。- 碰撞變得太大。以下代碼通過蠻力執行此操作，從小開始並將它加倍直到它足夠大。因此整個計算需要時間在哪裡是解決方案。人數的整個概率分佈被計算。# # 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) }舉個例子，人群中至少有 8 個人同一個生日更有可能是最少人數：, 由計算得出
birthday(7)。只需幾秒鐘。這是部分輸出的圖：
這個問題的一個特殊版本在將生日悖論擴展到超過 2 個人中得到解決，它涉及到 a 的情況。滾動很多次的雙面骰子。
代碼
# 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) }答案是通過
print(tmultinom(100,20,6), digits=15)0.267747907805267
