Simulating Dice

Simulating Dice#

Simulating random variables enables us to model uncertainty and variability in real-world systems and allows us to investigate randomness. In this section we will start simulating random variables in Python.

Tip

Activate Python on this page with the rocket icon () above, and experiment with changing the code.

import numpy as np

# If you want to be able to generate the same numbers every time, you can set a seed below
# np.random.seed(42)

The most basic building blocks are uniform random variables. We will gloss over the mechanics of how random number generators are produced and simply call them, using the numpy package. Below is our first simulated \(X\sim U(0,1)\) random variable.

X = np.random.uniform()

X

The function np.random.uniform can be supplied two parameters, a lower limit and an upper limit, which are by default set to \(0\) and \(1\). Above you can enter these parameter to produce instead a uniformly distributed real number between for instance \(-1\) and \(4\) by entering \(-1, 4\) in the brackets.

We can use this uniformly distributed random variable to simulate various other distributions. For instance we may simulate a six sided die:

X = int(np.random.uniform(1,7))

X

Indeed in the code above we randomly generate a uniformly distributed number between 1 and 7 then round down to the nearest integer, which produces a(n ideal) die. You may want to try to change the code to simulate a \(4\) or \(8\) sided die.

Very often we will need to simulate not one die, but many dice at the same time. The third parameter of the uniform distribution is the size. We may use this parameter to produce an array of many random numbers at once. Each random number is independently uniformly distributed.

X = np.random.uniform(1,7,10).astype(int)

X

As a check we may consider the relative frequency of each number when we throw many six sided dice. In the following we first define rolling some number of dice as a function for future use. Then we check whether the relative frequency is \(\frac{1}{6}=0.1666\ldots\) for each outcome.

# Simulates throwing a specified number of dice with a specified number of sides returning the outcomes as integers.
def throw_dice(n_sides = 6, n_dice = 1):
    return np.random.uniform(1, n_sides + 1, n_dice).astype(int)

n_dice = 100000
n_sides = 6
dice_throws = throw_dice(n_sides, n_dice)

for i in range(1, n_sides + 1): 
    print(f"Proportion of times we rolled a {i}: {np.count_nonzero(dice_throws==i)/n_dice:.3f}")

Finally we can now simulate not just individual die rolls, but also sums of dice rolls leading to simulations of non-uniformly distributed random variables. For example below we simulate many rolls of a pair of \(4\) sided die rolls and check that we obtain the expected distribution of sums.

# Simulate rounds of dice throws and return an array with the sum of outcomes for each round.
def throw_and_sum_dice(n_sides = 6, n_throws = 1, n_dice_per_throw = 2):
    return np.sum(throw_dice(n_sides, (n_throws, n_dice_per_throw)), axis = 1)

n_sides = 4
n_throws = 100000
n_dice_per_throw = 2
summed_dice_throws = throw_and_sum_dice(n_sides, n_throws, n_dice_per_throw)

for i in range(n_dice_per_throw, n_dice_per_throw * n_sides + 1):
    print(f"Proportion of times we rolled a {i}: {np.count_nonzero(summed_dice_throws==i)/n_throws:.3f}")

Bonus: Visualisation#

Here is a way to visualise the distributions of the (sums of) dice rolls.

import micropip
await micropip.install("seaborn")

import seaborn as sns
import matplotlib.pyplot as plt

sns.set_theme(style="white")
fig, axes = plt.subplots(1, 2, figsize=(10, 4))  # 1 row, 2 columns

sns.histplot(dice_throws, discrete=True, stat="density", edgecolor="black", ax=axes[0])
axes[0].set(xlabel="Outcome", ylabel="Frequency", title=f"Dice Throw Outcome Frequencies")

sns.histplot(summed_dice_throws, discrete=True, stat="density", edgecolor="black", ax=axes[1])
axes[1].set(xlabel="Summed Dice Outcome", ylabel="Frequency", title=f"Summed Dice Throws Frequencies \nfor {n_sides}-sided dice, with {n_dice_per_throw} dice per throw, for {n_throws} throws")

plt.tight_layout()  # Adjust layout to prevent overlapping
plt.show()
Contents