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

发布于 2024-07-09 07:11:31 字数 2221 浏览 8 评论 0原文

好吧,这个问题困扰了我好几年了。 如果你在学校学过统计学和高等数学,那么现在就转身吧。 为时已晚。

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

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

如果您编写一些代码(见下文),将这些骰子掷几百万次,然后计算收到每个数字作为最终结果的频率,您会得到一条在 1 左侧非常平坦的曲线,在 1 和 1 之间大约 45° 度。 60 且在 60 之上持平。掷出 30.5 或更好的机会大于 50%,掷出好于 18 的机会为 80%,掷出好于 0 的机会为 97%。

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

背景:对于我们的角色扮演游戏“星之丛林”,我们寻找一种方法来控制随机事件。 上面的规则保证了你尝试的东西会得到更稳定的结果:)

对于周围的极客来说,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

Okay, this bugged me for several years, now. If you sucked in statistics and higher math at school, turn away, now. Too late.

Okay. Take a deep breath. Here are the rules. Take two thirty sided dice (yes, they do exist) and roll them simultaneously.

  • Add the two numbers
  • If both dice show <= 5 or >= 26, throw again and add the result to what you have
  • If one is <= 5 and the other >= 26, throw again and subtract the result from what
    you have
  • Repeat until either is > 5 and < 26!

If you write some code (see below), roll those dice a few million times and you count how often you receive each number as the final result, you get a curve that is pretty flat left of 1, around 45° degrees between 1 and 60 and flat above 60. The chance to roll 30.5 or better is greater than 50%, to roll better than 18 is 80% and to roll better than 0 is 97%.

Now the question: Is it possible to write a program to calculate the exact value f(x), i.e. the probability to roll a certain value?

Background: For our role playing game "Jungle of Stars" we looked for a way to keep random events in check. The rules above guarantee a much more stable outcome for something you try :)

For the geeks around, the code in 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

如果你对这篇内容有疑问,欢迎到本站社区发帖提问 参与讨论,获取更多帮助,或者扫码二维码加入 Web 技术交流群。

扫码二维码加入Web技术交流群

发布评论

需要 登录 才能够评论, 你可以免费 注册 一个本站的账号。

评论(4

挥剑断情 2024-07-16 07:11:31

在我能够理解它之前,我必须先重写你的代码:

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 的 字典

不管怎样,抛开编程风格的问题不谈:假设 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),那么您将再次投掷,因此此项为 F( 的表达式贡献了 (1/30)(1/30)(5+F(x-5)) X)。 所以你会得到一些长得令人发指的表达式,比如

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

它是 900 多个项的总和,每个项对应 [30]×[30] 中的每一对 (a,b)。 两个 ≤ 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)]

形式的关系

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

(例如),您(从 F 和 G 的这些表达式)可以得到适当矩阵 A 和适当向量 B (可以计算)的 ,因此对于足够小的 x,从 V(x) = [0 0] 形式的初始值开始,您可以在您想要任意接近精度的范围内找到 x 的 F(x) 和 G(x)。 (而你的 f(x),即抛出 x 的概率,就是 F(x)-F(x-1),所以结果也出来了。)

可能有更好的方法。 话虽如此,但你为什么要这样做呢? 无论您想要哪种分布,都有很好且简单的概率分布,具有适当的参数,具有良好的特性(例如小方差、单边误差等)。 没有理由编写自己的临时程序来生成随机数。

I had to first rewrite your code before I could understand it:

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

Maybe you might find this less readable; I don't know. (Do check if it is equivalent to what you had in mind.) Also, regarding the way you use "result" in your code -- do you know of Python's dicts?

Anyway, matters of programming style aside: Suppose F(x) is the CDF of OW60(1), i.e.

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

Similarly let

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

Then you can calculate F(x) from the definition, by summing over all (30×30) possible values of the result of the first throw. For instance, if the first throw is (2,3) then you'll roll again, so this term contributes (1/30)(1/30)(5+F(x-5)) to the expression for F(x). So you'll get some obscenely long expression like

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

which is a sum over 900 terms, one for each pair (a,b) in [30]×[30]. The pairs (a,b) with both ≤ 5 or both ≥26 have a term a+b+F(x-a-b), the pairs with one ≤5 and one ≥26 have a term a+b+G(x-a-b), and the rest have a term like (a+b), because you don't throw again.

Similarly you have

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

Of course, you can collect coefficients; the only F terms that occur are from F(x-60) to F(x-52) and from F(x-10) to F(x-2) (for a,b≥26 or both≤5), and the only G terms that occur are from G(x-35) to G(x-27) (for one of a,b≥26 and the other ≤5), so there are fewer terms than 30 terms. In any case, defining the vector V(x) as

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

(say), you have (from those expressions for F and G) a relation of the form

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

for an appropriate matrix A and an appropriate vector B (which you can calculate), so starting from initial values of the form V(x) = [0 0] for x sufficiently small, you can find F(x) and G(x) for x in the range you want to arbitrarily close precision. (And your f(x), the probability of throwing x, is just F(x)-F(x-1), so that comes out as well.)

There might be a better way. All said and done, though, why are you doing this? Whatever kind of distribution you want, there are nice and simple probability distributions, with the appropriate parameters, that have good properties (e.g. small variance, one-sided errors, whatever). There is no reason to make up your own ad-hoc procedure to generate random numbers.

老旧海报 2024-07-16 07:11:31

我对 2000 万次投掷的样本做了一些基本统计。 结果如下:

Median: 17 (+18, -?) # This result is meaningless
Arithmetic Mean: 31.0 (±0.1)
Standard Deviation: 21 (+1, -2)
Root Mean Square: 35.4 (±0.7)
Mode: 36 (seemingly accurate)

误差是通过实验确定的。 算术平均值和众数非常准确,即使非常积极地更改参数似乎也不会对其产生太大影响。 我想中位数的行为已经得到了解释。

注意:不要将这些数字用于函数的正确数学描述。 使用它们可以快速了解分布情况。 对于其他任何事情,它们都不够准确(即使它们可能很精确。

也许这对某人有帮助。

编辑 2:

graph

仅基于 991 个值。 我本可以在其中塞入更多的价值观,但它们会扭曲结果。 这个样本恰好是相当典型的。

编辑 1:

以下是仅一个 60 面骰子的上述值,用于比较:

Median: 30.5
Arithmetic Mean: 30.5
Standard Deviation: 7.68114574787
Root Mean Square: 35.0737318611

请注意,这些值是计算得出的,而不是实验得出的。

I've done some basic statistics on a sample of 20 million throws. Here are the results:

Median: 17 (+18, -?) # This result is meaningless
Arithmetic Mean: 31.0 (±0.1)
Standard Deviation: 21 (+1, -2)
Root Mean Square: 35.4 (±0.7)
Mode: 36 (seemingly accurate)

The errors were determined experimentally. The arithmetic mean and the mode are really accurate, and changing the parameters even quite aggressively doesn't seem to influence them much. I suppose the behaviour of the median has already been explained.

Note: don't take these number for a proper mathematical description of the function. Use them to quickly get a picture of what the distribution looks like. For anything else, they aren't accurate enough (even though they might be precise.

Perhaps this is helpful to someone.

Edit 2:

graph

Based on just 991 values. I could've crammed more values into it, but they would've distorted the result. This sample happens to be fairly typical.

Edit 1:

here are the above values for just one sixty-sided die, for comparison:

Median: 30.5
Arithmetic Mean: 30.5
Standard Deviation: 7.68114574787
Root Mean Square: 35.0737318611

Note that these values are calculated, not experimental.

对你的占有欲 2024-07-16 07:11:31

复合无界概率是……不平凡的。 我本想以与 James Curran 相同的方式解决这个问题,但后来我从你的源代码中看到可能有第三组卷,第四组,依此类推。 这个问题是可以解决的,但远远超出了大多数骰子模拟器的范围。

是否有任何特殊原因需要从 -Inf 到 +Inf 的随机范围以及 1-60 左右如此复杂的曲线? 为什么2D30的钟形曲线不可接受? 如果您解释您的要求,很可能有人可以提供更简单且更有边界的算法。

Compound unbounded probability is... non-trivial. I was going to tackle the problem the same way as James Curran, but then I saw from your source code that there could be a third set of rolls, and a fourth, and so on. The problem is solvable, but far beyond most die rolling simulators.

Is there any particular reason that you need a random range from -Inf to +Inf with such a complex curve around 1-60? Why is the bell curve of 2D30 not acceptable? If you explain your requirements, it is likely someone could provide a simpler and more bounded algorithm.

您的好友蓝忘机已上羡 2024-07-16 07:11:31

好吧,走着瞧。 第二次投掷(有时会与第一次投掷相加或减去)在 31 左右有一个很好且易于预测的钟形曲线。当然,第一次投掷是问题所在。

对于第一卷,我们有 900 种可能的组合。

  • 50 种组合会导致添加第二卷。
  • 25 种组合会导致减去第二次掷骰。
  • 留下 825 个与第二卷钟形曲线相匹配的组合。

减法集(减法前)将在 (27..35) 范围内形成钟形曲线。
加法集的下半部分将在(2..10)范围内形成钟形曲线,而上半部分将在(52...60)范围内形成钟形曲线

我的概率有点生疏,所以我无法为您计算出确切的值,但应该清楚的是,这些值会导致可预测的值。

Well, let's see. The second throw (which will sometimes be added or subtracted to the first roll) has a nice easily predictable bell curve around 31. The first roll, of course, is the problem.

For the first roll, we have 900 possible combinations.

  • 50 combinations result in adding the second roll.
  • 25 combinations result in subtracting the second roll.
  • Leaving 825 combinations which match the bell curve of the second roll.

The subtracting set (pre-subtraction) will form a bell curve in the range (27..35).
The lower half of the adding set will form a bell curve in the range (2..10), while the upper half will form a bell curve in the range (52...60)

My probablity is a bit rusty, so I can't figure the exact values for you, but it should be clear that these lead to predictable values.

~没有更多了~
我们使用 Cookies 和其他技术来定制您的体验包括您的登录状态等。通过阅读我们的 隐私政策 了解更多相关信息。 单击 接受 或继续使用网站,即表示您同意使用 Cookies 和您的相关数据。
原文