Simulating Bertrand's Paradox

One of the most famous elementary problems in probability is Bertrand's Paradox. This stems from the seemingly simple question: what is the probability that the length of a random chord in a circle exceeds the length of the side of its inscribed equilateral triangle? To put the lengths into perspective, the side of an equilateral triangle inscribed in a circle of radius $r$ has length $\sqrt{3}r$. The YouTube channel Numberphile, in collaboration with Grant Sanderson from 3blue1brown, released an excellent video discussing this very problem.

Don't forget to watch the extra footage!

I won't dwell on the mathematical details here; the video does a great job of explaining the problem with gorgeous visuals. Instead, an enthusiastic programmer's first instinct screams "let's code it"; and python is more than up to the task. Putting the usual caveats of generating (pseudo)randomness using computers aside, the numpy.random module is a good place to start. With help from Rohan Kumar, here's some code illustrating three methods of generating random chords.

import numpy as np
from numpy import pi
from numpy.random import uniform


def uniform_angles(N):
    """Sample N chords by choosing two points uniformly from the circumference."""
    theta = uniform(0, 2 * pi, (2, N))
    delta = np.abs(theta[0] - theta[1])
    hits = np.logical_and(delta > 2 * pi / 3, delta < 4 * pi / 3)

    return hits


def uniform_polar(N):
    """Sample N chords by choosing the center in polar coordinates, uniform r and theta."""
    theta = uniform(0, 2 * pi, N)  # Redundant
    radius = uniform(0, 1, N)
    hits = radius < 0.5

    return hits


def uniform_square(N):
    """Sample at most N chords by choosing the center in Cartesian coordinates, uniform x and y.
    Points outside the circle are discarded.
    """
    coords = uniform(-1, 1, (2, N))
    within = coords[0] ** 2 + coords[1] ** 2 <= 1
    x, y = coords[0][within], coords[1][within]
    hits = x ** 2 + y ** 2 < 0.5 ** 2

    return hits


def main():
    for method in [uniform_square, uniform_angles, uniform_polar]:
        # Sample up to 1 million chords per method
        hits = method(1000000)
        count = np.sum(hits)
        total = len(hits)
        print(f"{count:10d} / {total:10d} = {count/total:.6f} ({method.__name__})")


if __name__ == "__main__":
    main()

A typical run gives an output like so. The fractions are the proportions of long chords among all chords.

	196461 /     785231 = 0.250195 (uniform_square)
    	333021 /    1000000 = 0.333021 (uniform_angles)
    	500594 /    1000000 = 0.500594 (uniform_polar)

The mathematics may be simple, but it's always reassuring to see it play out! To round things off, here are $1500$ chords generated using the first method (choose two points uniformly from the circumference and join them). The red chords are the ones which are longer than $\sqrt{3}r$. The ring of blue chords clearly shows that these short chords must avoid the inner circle of radius $r / 2$.