具有大 n (>100) 的骰子数学

Question

我保证这不仅仅是另一个掷骰子的家庭作业问题。我实现了一个函数来计算在掷骰子s时获得小于总和的概率n m。我的函数适用于 的小值，n但我发现大值的n. 见附图。任何人都了解发生了什么？

我的概率函数

从此数学堆栈交换中实现

probability <- function(s, m, n) {

  i <- 0:((s-1-n) / m)
  m^(-n) * sum((-1)^i * choose(n, i) * choose(s - 1 - i * m, n))

}

开始打破 ~ n > 80

n <- 1:90 # number of dice
m <- 6 # number of sides
s <- floor(mean(1:m)*n) # sum of faces
p <- mapply(probability, s = s, m = m, n = n)
plot(n, p, main = paste("probability of rolling less than floor(", mean(1:m),"* n) as sum of n rolls of a", m, "sided dice"))

Answer 1

正如对原始问题的评论中提到的，问题在于概率函数要求 R 计算非常大的数字 ( choose(80,40) = 1.075072e+23)，而我们正在达到 R 的数值精度限制。

另一种不涉及大量数字而是使用大量数字的替代方法是运行蒙特卡罗模拟。这会生成骰子掷总和的分布，并将观察到的总和与分布进行比较。运行时间会更长，但更容易做到，并且不会出现数值精度问题。

mc <- Vectorize(function(s, m, n, reps = 10000) {
  x <- replicate(reps, sum(sample(m, n, replace = TRUE)))
  ecdf(x)(s-1)
})



n <- 1:90 # number of dice
m <- 6 # number of sides
s <- floor(mean(1:m)*n) # sum of faces
analytic_prob <- mapply(probability, s = s, m = m, n = n)
mc_prob <- mapply(mc, s = s, m = m, n = n)


plot(n, analytic_prob, main = paste("probability of rolling less than floor(", mean(1:m),"* n) as sum of n rolls of a", m, "sided dice"),
     sub = "monte carlo in red")
points(n, mc_prob, col = "red")