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

Question

我正试图编写一些Python代码来计算n个滚动的m边公平骰子之和成为某个值或更大值的概率。输入是所有骰子卷的总和、所做的卷数和骰子上的边数。输出是掷骰子总数达到或超过该值的概率百分比。

我的计算基于我在这篇论文中找到的方程，并对任何被轧制多次的单面模具进行缩放：https://digitalscholarship.unlv.edu/cgi/viewcontent.cgi?article=1025

我做了一些“有效”的代码，但是当骰子有很多面时速度非常慢，所以只对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面的模具，你将不得不迭代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%

在线试试吧！