计算两个D30复投的精确结果

Question

好的,现在这让我困扰了好几年了.如果你在学校学习统计数据和提高数学水平,现在就转身离开.太晚了.

好的.深吸一口气.这是规则.取两个三十面骰子(是的,它们确实存在)并同时滚动它们.

• 添加两个数字
• 如果两个骰子显示<= 5或> = 26,则再次抛出并将结果添加到您拥有的内容中
• 如果一个<= 5且另一个> = 26,则再次抛出并从您拥有的结果中减去结果
• 重复直到> 5和<26!

如果你写了一些代码(见下文),将这些骰子滚动几百万次并计算你接收每个数字作为最终结果的频率,你会得到一个非常平坦的曲线,1和1之间约为45°.滚动30.5或更好的机会大于50%,滚动优于18是80%并且滚动优于0是97%.

现在的问题是:是否可以编写一个程序来计算的确切值f(x)的,即滚某个值的概率是多少？

背景:对于我们的角色扮演游戏"Jungle of Stars",我们寻找一种方法来控制随机事件.对于你尝试的东西,上面的规则保证了更稳定的结果:)

对于极客来说,Python中的代码:

import random
import sys

def OW60 ():
    """Do an open throw with a "60" sided dice"""
    val = 0
    sign = 1

    while 1:
        r1 = random.randint (1, 30)
        r2 = random.randint (1, 30)

        #print r1,r2
        val = val + sign * (r1 + r2)
        islow = 0
        ishigh = 0
        if r1 <= 5:
            islow += 1
        elif r1 >= 26:
            ishigh += 1
        if r2 <= 5:
            islow += 1
        elif r2 >= 26:
            ishigh += 1

        if islow == 2 or ishigh == 2:
            sign = 1
        elif islow == 1 and ishigh == 1:
            sign = -1
        else:
            break

        #print sign

    #print val
    return val

result = [0] * 2000
N = 100000
for i in range(N):
    r = OW60()
    x = r+1000
    if x < 0:
        print "Too low:",r
    if i % 1000 == 0:
        sys.stderr.write('%d\n' % i)
    result[x] += 1

i = 0
while result[i] == 0:
    i += 1

j = len(result) - 1
while result[j] == 0:
    j -= 1

pSum = 0
# Lower Probability: The probability to throw this or less
# Higher Probability: The probability to throw this or higher
print "Result;Absolut Count;Probability;Lower Probability;Rel. Lower Probability;Higher Probability;Rel. Higher Probability;"
while i <= j:
    pSum += result[i]
    print '%d;%d;%.10f;%d;%.10f;%d;%.10f' % (i-1000, result[i], (float(result[i])/N), pSum, (float(pSum)/N), N-pSum, (float(N-pSum)/N))
    i += 1

Answer 1

在我理解之前,我必须先重写代码:

def OW60(sign=1):
    r1 = random.randint (1, 30)
    r2 = random.randint (1, 30)
    val = sign * (r1 + r2)

    islow  = (r1<=5)  + (r2<=5)
    ishigh = (r1>=26) + (r2>=26)

    if islow == 2 or ishigh == 2:
        return val + OW60(1)
    elif islow == 1 and ishigh == 1:
        return val + OW60(-1)
    else:
        return val

也许你会发现这个可读性较差; 我不知道.(请检查它是否与您的想法相同.)另外,关于在代码中使用"结果"的方式 - 你知道Python的dict吗？

无论如何,编程风格的问题除外:假设F(x)是OW60(1)的CDF,即

F(x) = the probability that OW60(1) returns a value ? x.

同样地

G(x) = the probability that OW60(-1) returns a value ? x.

然后,您可以通过对第一次投掷结果的所有(30×30)个可能值求和来计算定义中的F(x).例如,如果第一次投掷是(2,3)那么你将再次滚动,所以这个术语贡献(1/30)(1/30)(5 + F(x-5))到F的表达式( X).所以你会得到一些如此猥亵的长期表达

F(x) = (1/900)(2+F(x-2) + 3+F(x-3) + ... + 59+F(x-59) + 60+F(x-60))

这是一个超过900项的总和,每对(a,b)一个[30]×[30].两个≤5或两者≥26的对(a,b)具有项a + b + F(xab),一个≤5且一个≥26的对具有项a + b + G(xab),并且其余的有一个像(a + b)这样的术语,因为你不再扔了.

你也一样

G(x) = (1/900)(-2+F(x-2) + (-3)+F(x-3) + ... + (-59)+F(x-59) + (-60)+F(x-60))

当然,你可以收集系数; 唯一出现的F项是从F(x-60)到F(x-52)和从F(x-10)到F(x-2)(对于a,b≥26或两者都是≤5),唯一出现的G项是从G(x-35)到G(x-27)(a中的一个,b≥26,另一个≤5),所以项数少于30项.在任何情况下,将向量V(x)定义为

V(x) = [F(x-60) G(x-60) ... F(x-2) G(x-2) F(x-1) G(x-1) F(x) G(x)]

(比如说),你有(从F和G的表达式)形式的关系

V(x) = A*V(x-1) + B

对于适当的矩阵A和适当的向量B(你可以计算),所以从形式V(x)= [0 0]的初始值开始,对于足够小的x,你可以找到F(x)和G(x )对于你想要任意关闭精度的范围内的x.(而你的f(x),投掷x的概率,只是F(x)-F(x-1),因此也是如此.)

可能有更好的方法.不过,所有人都说完了,你为什么要这样做？无论你想要什么样的分布,都有很好的和简单的概率分布,具有适当的参数,具有良好的性质(例如小方差,单侧误差,等等).没有理由构建自己的临时程序来生成随机数.