用 S 面掷 N 个骰子时 T 总眼睛的概率

Question

我想计算 n 骰子的所有眼睛与 s 边（编号从 1 到 s）的总和等于 t 的事件的概率。我的语言是 Python 3。

我目前的方法几乎是一种尝试计数的解决方案，并且只适用于少量（运行 probability(10, 10, 50) 已经吃掉了我所有的 RAM 并迫使我硬重置）：

import itertools
def probability(n, s, t):
    all_rolls=list(itertools.product(range(1,s+1), repeat=n))
    target_rolls=[l for l in all_rolls if sum(l)==t]
    return round(len(target_rolls)/len(all_rolls), 4)

但老实说，我不知道如何解决这个问题。你能帮我走上正轨吗？

Answer 1

itertools.product 对于大量边数 > 5 和边数 > 6 来说太慢了。在我的机器上有 dice_number: 10 和边数: 10 需要一个半小时来计算。 相反，我使用numpy.polypow 函数来计算目标，并且计算时间不到一秒。

from numpy.polynomial.polynomial import polypow

def probability(dice_number, sides, target):
    """
    Using numpy polynomial
    The number of ways to obtain x as a sum of n s-sided dice
    is given by the coefficients of the polynomial:

    f(x) = (x + x^2 + ... + x^s)^n
    """

    # power series (note that the power series starts from x^1, therefore
    # the first coefficient is zero)
    powers = [0] + [1] * sides
    # f(x) polynomial, computed used polypow in numpy
    poly = polypow(powers, dice_number)
    return poly[target] / sides ** dice_number if target < len(poly) else 0

Answer 2

首先：可能的滚动组合总数将始终为s**n，因此您无需存储所有可能性的列表即可获得其长度。同样，您可以只保留所需结果的总数，而不是保留它们的列表以节省内存空间，但它仍然不会大大加快函数的速度：

def probability(n, s, t):
    all_rolls = itertools.product(range(1,s+1), repeat=n) #no list, leave it a generator
    target_rolls=sum(1 for l in all_rolls if sum(l)==t) #just total them up
    return round(target_rolls/s**n, 4)

一种更有效的计算可能性的方法是使用dict 和一些巧妙的迭代。每个字典将使用滚动值作为键和频率作为值，每次迭代 prev 将是前一个 X 骰子的此字典，cur 将通过添加另一个骰子从中更新：

import collections
def probability(n, s, t):
    prev = {0:1} #previous roll is 0 for first time
    for _ in range(n):
        cur = collections.defaultdict(int) #current probability
        for r,times in prev.items():
            for i in range(1,s+1):
                #if r occured `times` times in the last iteration then
                #r+i have `times` more possibilities for the current iteration.
                cur[r+i]+=times
        prev = cur #use this for the next iteration

    return cur[t] / s**n
    #return round(cur[t] / s**n , 4)

注意 1：因为cur 是一个默认字典，试图查找给定输入无法查找的数字将返回 0

注意 2：由于此方法将包含所有可能结果的字典放在一起，您可以返回 cur 并在同一次掷骰上计算多个不同的可能结果。

Answer 3

停止列清单。只需使用惰性求值即可。

from itertools import product

def prob(dice, pips, target):
    rolls = product(range(1, pips+1), repeat=dice)
    targets = sum(1 for roll in rolls if sum(roll) == target)
    return targets / pips**dice

测试：

for i in range(5, 26):
    print(i, prob(5, 5, i))
print('sum: ', sum(prob(5, 5, i) for i in range(5, 26)))
# prints
5 0.00032
6 0.0016
7 0.0048
8 0.0112
9 0.0224
10 0.03872
11 0.0592
12 0.0816
13 0.1024
14 0.1168
15 0.12192  # symmetrical, with peak in middle
16 0.1168
17 0.1024
18 0.0816
19 0.0592
20 0.03872
21 0.0224
22 0.0112
23 0.0048
24 0.0016
25 0.00032
sum:  1.0000000000000002

编辑：删除未使用的定义