dice-rolling-simulation.py

import numpy as np

n_rolls = 10000000
n_sides = 20

rolls_1 = np.random.randint(n_sides, size=(n_rolls, 2)) + 1
rolls_2 = np.random.randint(n_sides, size=(n_rolls, 2)) + 1

print('-' * 30)
print('advantage of disadvantages')
print("expected:  ", (2 * n_sides ** 2 - 1) / n_sides ** 4)
print("simulation:", np.mean(np.maximum(rolls_1.min(axis=1), rolls_2.min(axis=1))
 == n_sides))

print('-' * 30)
print('disadvantage of advantages')
print("expected:  ", ((2 * n_sides - 1) / n_sides ** 2) ** 2)
print("simulation:", np.mean(np.minimum(rolls_1.max(axis=1), rolls_2.max(axis=1))
 == n_sides))

print('-' * 30)
print('single roll')
print("expected:  ", 1 / n_sides)
print("simulation:", np.mean(rolls_1[:, 0] == n_sides))