如何计算n个滚动的m面骰子之和为某个值或更大的概率

Question

我正在尝试编写一些 Python 代码，该代码将计算 n 掷出的 m 面公平骰子之和为某个值或更大的概率。输入是所有掷骰子的总和、掷骰子的数量和骰子的边数。输出是骰子总和等于或大于该值的百分比机会。

我的计算基于我在本文中找到的方程式，并针对任何滚动任意次数的双面骰子进行缩放：https://digitalscholarship.unlv.edu/cgi/viewcontent.cgi?article=1025&context=grrj

我编写了一些“有效”的代码，但是当骰子有很多面时速度非常慢，所以只对 20 面或更少面的骰子有用。

import numpy as np

def probability_calculator(roll_total, num_of_rolls, dice_faces):       

    if num_of_rolls == 0:
        probability = 100    
    elif num_of_rolls == 1:
        probability = (dice_faces + 1 - roll_total) * 100/dice_faces    
    else:
        inverse = (dice_faces ** num_of_rolls) ** -1
        side_list = np.linspace(1, dice_faces, dice_faces)
        expanded_list = np.zeros(dice_faces * num_of_rolls)         
        stacked_side_list = side_list   

        for i in range(num_of_rolls - 1):
            stacked_side_list = np.vstack((stacked_side_list, side_list))
            
        index_array = np.zeros(num_of_rolls, dtype=int)
            
        while True:                
            value = 0

            for i in range(num_of_rolls):                    
                value = value + stacked_side_list[i][index_array[i]]
                
           expanded_list[int(value) - 1] += 1

           if sum(index_array) == (dice_faces - 1) * num_of_rolls:
               break

           for i in range(num_of_rolls):  
               if index_array[i] == dice_faces - 1:
                   index_array[i] = 0          
               else:
                    index_array[i] += 1                        
                    break

       probability = inverse * sum(expanded_list[roll_total - 1:]) * 100
  
   return probability

如您所见，这是非常低效的代码，如果您只滚动四个 100 面骰子，您将不得不遍历 while 循环 100^4 = 100,000,000 次.....

我很确定有一些数学方程式可以简化此代码并使其运行速度提高许多数量级，但数学不是我最喜欢的科目，我不知道有任何方程式或 Python 函数可以提供帮助。

Answer 1

快速浏览了论文并被公式吓到了。实施了一种可能不同的方式，只需计算达到不同总数的频率：

from functools import lru_cache

@lru_cache(None)
def sum_freq(total, rolls, faces):
        if not rolls:
            return not total
        return sum(sum_freq(total - die, rolls - 1, faces)
                   for die in range(1, faces + 1))

def probability_calculator(roll_total, num_of_rolls, dice_faces):
    return sum_freq(roll_total, num_of_rolls, dice_faces) / dice_faces**num_of_rolls

“四个 100 面骰子”的演示：

prob_314 = probability_calculator(314, 4, 100)
prob_any = sum(probability_calculator(total, 4, 100)
               for total in range(1, 401))
print(f'{prob_314:%}')
print(f'{prob_any:%}')

输出：

0.113564%
100.000000%

10 100 面骰子的输出：

0.050065%
100.000000%

42 个这样的 100 面骰子的输出：

0.000000%
100.000000%

Try it online!