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

Question

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

我的概率函数

从此实现 math stack exchange

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")