Lesson 3: Randomness and counting

"""Galton board: every ball goes left or right at random at each of 12 rows of pegs.
Run: python galton-board.py   (needs: pip install numpy matplotlib)
"""
from math import comb

import numpy as np
import matplotlib.pyplot as plt

plt.style.use("dark_background")

rows, balls, p = 12, 5000, 0.5  # p: the chance of bouncing right
rng = np.random.default_rng()

bins = (rng.random((balls, rows)) < p).sum(axis=1)  # a ball's bin = how often it went right
counts = np.bincount(bins, minlength=rows + 1)

k = np.arange(rows + 1)
expected = [balls * comb(rows, i) * p**i * (1 - p) ** (rows - i) for i in k]  # the binomial curve

plt.figure(figsize=(7, 4))
plt.bar(k, counts, color="#5eead4", label="balls")
plt.xlabel("bin = times the ball went right")
plt.plot(k, expected, color="white", label="C(n, k) p^k (1-p)^(n-k)")
plt.legend()
plt.show()

# Try: p = 0.65 tilts the board; rows = 24 makes the bell narrower for its width

No Python on this computer? The notebook runs the same script, the exercises and the answers in your browser, on Google Colab.

Line by line

1. rng = np.random.default_rng()

   numpy’s random number generator. Give it a number, default_rng(1), and every run gives the same “random” numbers.

2. bins = (rng.random((balls, rows)) < p).sum(axis=1)

   A table of 5,000 × 12 random numbers between 0 and 1. Each one below p is True, a bounce to the right; adding up a row counts its Trues, which is the ball’s bin.

3. counts = np.bincount(bins, minlength=rows + 1)

   How many balls landed in each bin from 0 to 12; minlength makes sure the empty edge bins are counted as 0.

4. expected = [balls * comb(rows, i) * p**i * (1 - p) ** (rows - i) for i in k]

   A list comprehension: the binomial formula for every bin i, times the number of balls.

5. plt.bar(k, counts, color="#5eead4", label="balls")

   One bar per bin. The white line drawn after it is the formula, to compare.

Try it yourself

1. Run the script twice. Are the bars the same? Then make every run the same.

   Hint

   Give default_rng a number.

   Answer

   The bars differ a little each time. rng = np.random.default_rng(1) gives the same run every time.

2. How many of 5,000 balls should land in the centre bin of 12 rows?

   Hint

   The share of the centre bin is comb(12, 6) / 2**12.

   Answer

   5000 × 924 ÷ 4096 ≈ 1128.

3. Tilt the board with p = 0.65. Where is the bell’s centre now?

   Hint

   The centre is rows · p.

   Answer

   12 × 0.65 = 7.8, so the peak moves to bins 7 and 8.

Quiz