大数值（>100）的骰子点数统计数学问题

Question

我保证这不仅仅是又一个骰子点数作业问题。我实现了一个函数来计算掷出n枚m面骰子，总和小于值s的概率。我的函数可以计算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"))

- Lief Esbenshade

---

雪上加霜的是，我的朋友已经在Mathematica中实现了相同的算法，并且在处理大量n值时没有遇到任何问题。 - Lief Esbenshade

那些“选择”数字变得非常大。例如choose(80,40)。你的公式在数值上不稳定。也许尝试在对数尺度上计算会更好。 - MrFlick

当n很大时，choose会完全失去精度。也许你可以阅读 https://dev59.com/r1kR5IYBdhLWcg3w_CD-#40527881 以获取替代方案。 - ThomasIsCoding

谢谢，我正在尝试实现拉马努金逼近算法，但是在向量化函数和让它使用base::choose处理小值的n和k时遇到了问题。我会编辑问题以包含我的进展。 - Lief Esbenshade

已解决NaN问题，但是当n > 80时精度仍然会出现问题。 - Lief Esbenshade

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

- Lief Esbenshade

Answer 2

问题是由于R的数值精度限制引起的。正如评论者所指出的那样，我计算的n个k值实在太大了（choose(80,40) = 1.075072e+23）。

我们可以使用对数来尝试将问题保持在R的计算限制内。这是拉马努金方法的实现。不幸的是，近似误差会累积，并且精度下降得更快。概率函数需要添加和减去一系列非常大的数字，以获得介于0和1之间的最终值，并且不能容忍任何不精确性。

0）重写概率函数以分步进行

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

  # Probability of getting less than s
  i <- 0:((s-1-n) / m)

  c1 <- choose(n, i)
  c2 <- choose(s - 1 - i * m, n)

  seq <- (-1)^i * (c1 * c2)

  m^(-n) * sum(seq)

}

1) 实现log(x!)的近似值

# using the 'ramanujan' method
ramanujan <- function(n){
  n * log(n) - n + log(n * (1 + 4*n * (1 + 2*n))) / 6 + log(pi) / 2
}

# confirm Ramanujan works correctly
n <- 1:200
diff <- log(factorial(n)) - ramanujan(n)
plot(n, diff) # r returns inf for factorial(171), but up to there the numbers match

2) 使用对数逼近重写choose函数。

#' This function returns log(choose(n,k)) 
log_nck <- Vectorize(function(n, k) {
  if(n <= k | n < 1 | k < 1) return(log(choose(n,k))) # logs don't like 0 or neg numbers

  return((ramanujan(n) - ramanujan(k) - ramanujan(n-k)))
})

# Check that choose function works
n <- seq(10, 100, 10)
k <- seq(5, 50, 5)
c_real <- log(choose(n, k))
c_approx <- log_nck(n, k)
# If we print them, they appear to match
print(c_real)
print(c_approx)
# and the difference shows pretty small errors. 
print(c_real - c_approx)

3) 使用对数组合重写概率函数。

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

  # Probability of getting less than s
  i <- 0:((s-1-n) / m)

  c1 <- log_nck(n, i)
  c2 <- log_nck(s - 1 - i * m, n)

  seq <- (-1)^i * exp(c1 + c2)

  return(m^(-n) * sum(seq))

}

最终测试

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)
newp <- mapply(new_probability, s = s, m = m, n = n)

plot(n, p, main = "Original in black, approximation in red")
points(n, newp, col = "red")

- Lief Esbenshade