python - python中的蒙特卡洛模拟骰子

Question

我有一个与基于掷 2 个骰子的概率的蒙特卡罗模拟相关的问题。当用 python 编码时，如何呈现总和大于 n 且小于 m 的事实？作为一个例子，我在 mathlab 中做了这个:

NT = 10^5; %number of throws
log = zeros(1,12);
for throw = 1:NT
    dices = ceil(6*rand(1,2));
    s = sum(dices);
    log(s) = log(s)+1;
end
p = 36*log(6:9)/NT;
s1 = sum(round(p))

在上面的示例中，我假设 n 为 5，m 为 10。

谢谢

Answer 1

在每个循环中，您想要模拟两次单独的随机骰子 throw 。我的以下代码片段使用 list (如果愿意，您可以使用dict)来存储NT的结果模拟:

import random

num_throws = 10**5  # NT
roll_log = [0] * 12  # Generate list for dice roll tallies

for i in range(num_throws):
    # Random integer between 1 and 6 inclusive for each dice
    dice_1 = random.randint(1, 6)
    dice_2 = random.randint(1, 6)

    # Sum the random dice and increment the tally for that particular roll total
    roll_sum = dice_1 + dice_2
    roll_log[roll_sum-1] += 1  # minus 1 because Python is 0-indexed

要处理结果数据，您可以通过roll_log[roll-1]访问结果列表中特定掷骰子的计数。哪里2 <= roll <= 12 (roll = 1 的概率为零，因为 2 个以上的骰子是不可能的)。以下 for 循环只是如何访问 NT 结果的示例。如果您不熟悉 Python 中的枚举，请进行模拟:

for i, tally in enumerate(roll_log):
    roll_prob = float(tally) / num_throws  # Experimental probability of roll
    roll = i + 1  # Since Python lists are 0-indexed
    print('{}: {}/{} = {}'.format(roll, tally, num_throws, roll_prob))

输出:

1: 0 / 100000 = 0
2: 2741 / 100000 = 0.02741
3: 5518 / 100000 = 0.05518
4: 8202 / 100000 = 0.08202
5: 11235 / 100000 = 0.11235
6: 14046 / 100000 = 0.14046
7: 16520 / 100000 = 0.1652
8: 13799 / 100000 = 0.13799
9: 11025 / 100000 = 0.11025
10: 8459 / 100000 = 0.08459
11: 5672 / 100000 = 0.05672
12: 2783 / 100000 = 0.02783

专门解决您问题的最后一部分，找到掷骰子在 n = 5 之间的概率。和m = 10不包含在内，这可以使用称为列表切片的方法来完成:

n = 5
m = 10
                                 #    6      7      8      9
rolls_between = roll_log[n:m-1]  # [14046, 16520, 13799, 11025]
sum_rolls_between = sum(rolls_between)  # 55390
prob_between = float(sum_rolls_between) / num_throws  # 0.5539

注意:float sum_rolls_between 的转换或num_throws最后一行对于获得十进制输出至关重要，因为在应用数学 floor() 后，Python 中两个整数之间的除法总是会产生整数输出。功能。换句话说，如果不将这两个值之一更改为浮点值，结果将为 0。