掷 n 个骰子后获得特定总和的概率。Ruby

Question

寻找n个骰子的点数和概率，最佳解决方案是什么？

1. 求平均值。
2. 求标准差。
3. 对于小于x的数字计算z-score。
4. 对于大于x的数字计算z-score。
5. 将两者都转换为概率。
6. 用其中一个概率减去另一个概率。

这是我迄今为止所做的。

# sides - number of sides on one die
def get_mean(sides)
  (1..sides).inject(:+) / sides.to_f
end

def get_variance(sides)
  mean_of_squares = ((1..sides).inject {|sum, side| sum + side ** 2}) / sides.to_f
  square_mean = get_mean(sides) ** 2

  mean_of_squares - square_mean
end

def get_sigma(variance)
  variance ** 0.5
end

# x - the number of points in question
def get_z_score(x, mean, sigma)
  (x - mean) / sigma.to_f
end

# Converts z_score to probability
def z_to_probability(z)
  return 0 if z < -6.5
  return 1 if z > 6.5

  fact_k = 1
  sum = 0
  term = 1
  k = 0

  loop_stop = Math.exp(-23)
  while term.abs > loop_stop do
    term = 0.3989422804 * ((-1)**k) * (z**k) / (2*k+1) / (2**k) * (z**(k+1)) / fact_k
    sum += term
    k += 1
    fact_k *= k
  end

  sum += 0.5
  1 - sum
end

# Calculate probability of getting 'х' total points by rolling 'n' dice with 'sides' number of sides.
def probability_of_sum(x, n, sides=6)

  mean = n * get_mean(sides)
  variance = get_variance(sides)
  sigma = get_sigma(n * variance)

  # Rolling below the sum
  z1 = get_z_score(x, mean, sigma)
  prob_1 = z_to_probability(z1)

  # Rolling above the sum
  z2 = get_z_score(x+1, mean, sigma)
  prob_2 = z_to_probability(z2)

  prob_1 - prob_2
end

# Run probability for 100 dice
puts probability_of_sum(400, 100)

我关心的是，当我选择x = 200时，概率为0。这个解决方案正确吗？

- Pav31

---

如果我理解正确的话，不是的。如果你掷100次骰子，所有的点数都是2，总和将是200。因此，这种情况发生的概率是存在的。 - Mangesh Tamhankar

这也是我想的。然而，根据68-95-99.7（经验）法则，或3西格玛规则，“大约68%从正态分布中抽取的值在距离平均值1个标准差_σ_内; 大约95%的值在2个标准差内; 而大约99.7%在3个标准差内。”在这种情况下，使用100个6面骰子，平均值= 350和σ = 17。这意味着99.7%的值将落在299到401的范围内。（350 +/- 17 * 3） - Pav31

还有一点需要注意的是，该算法存在一些根本性问题。即使是运行像 x=12 和 n=2 这样非常简单的计算，其结果应该是 1/36，但实际上并不正确。我认为 @NeilSlater 关于离散性方面的看法可能是正确的。 - Mangesh Tamhankar

在像掷两个骰子这样简单的情况下，“蒙特卡罗”模拟是更好的解决方案。但这种情况不同。 - Pav31

蒙特卡罗方法并不是一个关注结果准确性的好解决方案。对于像相同骰子之和这样简单的分布，可以以完美精度计算，并且比任何合理误差的蒙特卡罗方法更快。 - Neil Slater

Answer 1

有一个精确的解决方案，涉及二项式系数的交替求和。我已经在一些地方写出来了（在Quora和MSE上），你也可以在其他地方找到它，尽管有一些有缺陷的版本。请注意，如果你实现它，你可能需要取消比最终结果大得多的二项式系数，如果你使用浮点运算，你可能会失去太多的精度。

Neil Slater建议使用动态规划计算卷积是一个好主意。它比二项式系数的求和慢，但相当健壮。你可以通过几种方式加速它，例如使用平方指数法，以及使用快速傅里叶变换，但许多人会觉得这些方法过于复杂。

为了修正你的方法，你应该使用正态近似中的（简单）连续性校正，并且限制在有足够的骰子且在远离最大值和最小值时进行评估的情况下，您期望正态近似在绝对或相对意义上是好的。连续性校正是将计数n替换为从n-1/2到n+1/2的区间。

总共掷出200点的精确次数为7745954278770349845682110174816333221135826585518841002760次，因此概率为除以6^100，大约为1.18563 x 10^-20。

带有简单连续性校正的正态近似为Phi((200.5-350)/sqrt(3500/12))-Phi((199.5-350)/sqrt(3500/12)) = 4.2 x 10^-19。这在绝对值方面是准确的，因为它非常接近于0，但由于偏差达到35倍，所以在相对值方面并不好。在靠近中心的地方，正态近似提供了更好的相对近似。

- Douglas Zare

Answer 2

将两个独立概率分布的结果相加，等同于将这两个分布进行卷积。如果这些分布是离散的，则称为离散卷积。

因此，如果单个骰子表示为：

probs_1d6 = Array.new(6) { Rational(1,6) }

然后可以这样计算2d6：

probs_2d6 = []
probs_1d6.each_with_index do |prob_a,i_a|  
  probs_1d6.each_with_index do |prob_b,i_b| 
    probs_2d6[i_a + i_b] = ( probs_2d6[i_a + i_b] || 0 ) + prob_a * prob_b
  end
end

probs_2d6
# =>  [(1/36), (1/18), (1/12), (1/9), (5/36), (1/6), 
#      (5/36), (1/9), (1/12), (1/18), (1/36)]

虽然这对于骰子的面数是n平方的，而且一个完全合理的组合可以减少这个数量，但这样做通常不太灵活，适用于更复杂的设置。这种方法的好处在于您可以继续添加更多的骰子并进行其他更奇特的组合。例如，要获得4d6，您可以将2d6的两个结果卷积在一起。使用有理数允许您避开浮点精度问题。

我跳过了一个细节，您确实需要存储初始偏移量（对于普通的六面骰子为+1），并将其加在一起，以便知道概率匹配到什么位置。

我在游戏骰子宝石中制作了一个更复杂的逻辑，使用浮点数而不是有理数，可以处理一些其他骰子组合。

以下是您的方法的基本重写，使用上述方法以天真的方式（仅逐个组合每个骰子的效果）：

def probability_of_sum(x, n, sides=6)
  return 0 if x < n
  single_die_probs = Array.new(sides) { Rational(1,sides) }

  combined_probs = [1] # Represents chance of getting 0 when rolling 0 dice :-)

  # This is not the most efficient way to do this, but easier to understand
  n.times do
    start_probs = combined_probs
    combined_probs = []
    start_probs.each_with_index do |prob_a,i_a|  
        single_die_probs.each_with_index do |prob_b,i_b| 
          combined_probs[i_a + i_b] = ( combined_probs[i_a + i_b] || 0 ) + prob_a * prob_b
        end
    end
  end

  combined_probs[ x - n ] || 0
end

puts probability_of_sum(400, 100).to_f
# => 0.0003172139126369326

请注意，该方法实际上计算了100到600的完整概率分布，因此您只需要调用它一次并将数组（加上偏移量+100）存储一次，就可以执行其他有用的操作，例如获取大于某个数字的概率。由于在Ruby中使用了数，所有这些操作都具有完美的精度。

因为在您的情况下，您只有一种类型的骰子，所以我们可以避免使用，直到最后，仅使用整数（基本上是组合值的计数），并将其除以组合总数（面数的次幂）。这样更快，并且在不到一秒钟内返回100个骰子的值。

def probability_of_sum(x, n, sides=6)
  return 0 if x < n
  combined_probs = [1] # Represents chance of getting 0 when rolling 0 dice :-)

  n.times do
    start_probs = combined_probs
    combined_probs = []
    start_probs.each_with_index do |prob_a,i_a|  
        sides.times do |i_b| 
          combined_probs[i_a + i_b] = ( combined_probs[i_a + i_b] || 0 ) + prob_a
        end
    end
  end

  Rational( combined_probs[ x - n ] || 0, sides ** n )
end

- Neil Slater

---

只是一个小澄清。独立变量之和的概率密度或质量函数与密度或质量函数的卷积不是类似的，而是确切的卷积。 - Robert Dodier

Answer 3

这是我的最终版本。

1. 在 get_z_score 中改变了总和的偏移量，分别为 x-0.5 和 x+0.5，以获得更精确的结果。
2. 添加了 return 0 if x < n || x > n * sides，以覆盖情况，其中总和小于骰子数量，大于骰子数乘以骰子面数。
3. 添加了基准测试及其结果。

主要功能

# sides - number of sides on one die
def get_mean(sides)
  (1..sides).inject(:+) / sides.to_f
end

def get_variance(sides)
  mean_of_squares = ((1..sides).inject {|sum, side| sum + side ** 2}) / sides.to_f
  square_mean = get_mean(sides) ** 2

  mean_of_squares - square_mean
end

def get_sigma(variance)
  variance ** 0.5
end

# x - the number of points in question
def get_z_score(x, mean, sigma)
  (x - mean) / sigma.to_f
end

# Converts z_score to probability
def z_to_probability(z)
  return 0 if z < -6.5
  return 1 if z > 6.5

  fact_k = 1
  sum = 0
  term = 1
  k = 0

  loop_stop = Math.exp(-23)
  while term.abs > loop_stop do
    term = 0.3989422804 * ((-1)**k) * (z**k) / (2*k+1) / (2**k) * (z**(k+1)) / fact_k
    sum += term
    k += 1
    fact_k *= k
  end

  sum += 0.5
  1 - sum
end

# Calculate probability of getting 'х' total points by rolling 'n' dice with 'sides' number of sides.
def probability_of_sum(x, n, sides=6)
  return 0 if x < n || x > n * sides

  mean = n * get_mean(sides)
  variance = get_variance(sides)
  sigma = get_sigma(n * variance)

  # Rolling below the sum
  z1 = get_z_score(x-0.5, mean, sigma)
  prob_1 = z_to_probability(z1)

  # Rolling above the sum
  z2 = get_z_score(x+0.5, mean, sigma)
  prob_2 = z_to_probability(z2)

  prob_1 - prob_2
end

基准测试

require 'benchmark'

Benchmark.bm do |x|
  x.report { @prob = probability_of_sum(350, 100).to_f }
  puts "\tWith x = 350 and n = 100:"
  puts "\tProbability: #{@prob}"
end
puts

Benchmark.bm do |x|
  x.report { @prob = probability_of_sum(400, 100).to_f }
  puts "\tWith x = 400 and n = 100:"
  puts "\tProbability: #{@prob}"
end
puts

Benchmark.bm do |x|
  x.report { @prob = probability_of_sum(1000, 300).to_f }
  puts "\tWith x = 1000 and n = 300:"
  puts "\tProbability: #{@prob}"
end

Result

       user     system      total        real
   0.000000   0.000000   0.000000 (  0.000049)
    With x = 350 and n = 100:
    Probability: 0.023356331366255034

       user     system      total        real
   0.000000   0.000000   0.000000 (  0.000049)
    With x = 400 and n = 100:
    Probability: 0.00032186531055478085

       user     system      total        real
   0.000000   0.000000   0.000000 (  0.000032)
    With x = 1000 and n = 300:
    Probability: 0.003232390001131513

- Pav31

---

嗨，感谢解决方案！你能解释一下0.3989422804的含义以及为什么loop_stop = Math.exp(-23)吗？ - AlexMrKlim

它取自这里，该问题源自于这里。您还可以尝试搜索“z分数到概率公式”以获得更好的理解。 - Pav31

Answer 4

我也使用蒙特卡罗方法解决了这个问题，结果相对接近。

# x - sum of points to find probability for
# n - number of dice
# trials - number of trials
def monte_carlo(x, n, trials=10000)
  pos = 0

  trials.times do
    sum = n.times.inject(0) { |sum| sum + rand(1..6) }
    pos += 1 if sum == x
  end

  pos / trials.to_f
end

puts monte_carlo(300, 100, 30000)

- Pav31