求 n 个骰子滚动总和的概率的最佳解决方案是什么? 我正在通过查找来解决它
x
x
这就是我到目前为止所做的。
# 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。
这是正确的解决方案吗?
将两个独立概率分布的结果相加与对两个分布进行卷积相同。如果分布是离散的,那么它是一个离散卷积。
因此,如果单个骰子表示为:
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)并将其加在一起,以便知道概率匹配什么。
我在gem games_dice中制作了这个逻辑的更复杂版本,采用浮点而不是Rational,可以处理其他一些骰子组合。
这是使用上述方法以一种简单的方式对您的方法进行基本重写(简单地一次组合一个骰子的效果):
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 中使用了
Rational
数字,所有这些都具有完美的精度。
因为在您的情况下,您只有一种类型的骰子,所以我们可以避免使用
Rational
直到最后,仅使用整数(本质上是组合值的计数),然后除以组合总数(幂的边数)卷数)。这要快得多,在一秒内返回 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
存在涉及二项式系数交替和的精确解。我已经在几个地方写了它(在 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 倍,所以从相对值来看,它并不是很好。正态近似给出了更接近中心的更好的相对近似值。
这是我的最终版本。
get_z_score
中的总和偏移分别更改为 x-0.5
和 x+0.5
,以获得更精确的结果。 return 0 if x < n || x > n * sides
以涵盖总和小于骰子数量且大于骰子数量乘以边数的情况。主要功能
# 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
结果
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
我也用Monte Carlo方法解决了这个问题,结果也比较接近。
# 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)
这是我长期以来一直在数学上问的一个问题。然而,我的数学不太好。最近我找到了这个公式。我不知道它是如何用代码编写的,但你可以在 Omnicalculator 网站
上找到它