Python implementation of a Galton board simulation

Question

I watched a YouTube show the other day in which a variant of a Galton board is used to make a random selection from a range of movies that will be reviewed in that show. In contrast to a standard Galton board in which the ball is usually dropped from the mid-point, which results in a normal distribution, the starting position of the ball in the show is random. While watching, I was wondering whether this change really resulted in a uniform distribution of choices, and I sat down to write a Galton board simulator to test.

I've ended up with a class Galton that can be initialized with different board parameters (number of rows, number of bins). It provides the method simulate() that can be called with a parameter controlling the starting position of each bead during a simulation run; the number of beads per simulation can also be specified.

I believe my implementation is correct, but apart from a visual inspection of the the histograms I have no good idea how to test that formally. The class and the methods contain doc-strings, and I've used type-hints (even if I'm not not a huge fan of this feature).

In which ways can my code be improved with regard to coding style, documentation style, or program logic? Is there something that I've missed?

from collections import Counter
from matplotlib import pyplot as plt
from matplotlib.ticker import PercentFormatter
from pandas import Series
from random import randint


class Galton:
    """
    A Galton board class that can be used to simulate the distribution into
    bins of a number of beads running through the board.

    The board is assumed to comprise an even number of 'row pairs'. Whenever a
    bead passes through a member of a row pair in a simulated run, the bin
    position of the bead is moved randomly either a half-bin unit to the
    left or to the right. Consequently, after passing a row pair, the bin
    position of the bead has either moved a full bin unit to the left, stayed
    at the same bin unit, or moved a full bin unit to the right. If this
    results in a bin position outside the board, the bead is bounced back into
    the board by one full bin unit.
    """
    def __init__(self, row_pairs: int, bins: int):
        """
        Initialize a Galton board with the specified number of row pairs and bins.

        Arguments
        --------
        row_pairs: int
            The number of "row pairs"
        bins: int
            The number of bins (zero is not included as a bin)
        """
        self.bins = bins
        self.rows = row_pairs * 2

    def is_valid(self, position: float) -> bool:
        """
        Check if the argument is a valid board position.

        Argument
        --------
        position: float
            The position value to be checked (in bin units, including
            half-bins)

        Returns
        -------
        check: bool
            True if the position value is valid, or false otherwise
        """
        return 0 < position < self.bins + 1

    def move_down(self, position: float) -> float:
        """
        Move the bead down one row on the board by updating its bin position.

        Argument
        --------
        position: float
            The current position of the bead (in bin units, including
            half-bins)

        Returns
        -------
        position: float
            The new position of the bead after passing a row (in bin units,
            including half-bins)
        """
        d = -0.5 if randint(0, 1) else 0.5
        position += d
        if not self.is_valid(position):
            position -= d * 2
        return position

    def run_bead(self, start=None) -> int:
        """
        Run a bead from its starting position into a final bin by passing it
        through the board.

        Argument
        --------
        start:
            The starting position of the bead (see `simulate()` for details)

        Returns
        -------
        position: int
            The bin in which the bead ends up after running through the board
        """
        position = start or randint(1, self.bins)
        if not self.is_valid(position):
            raise ValueError("Bin position out of range")
        for _ in range(self.rows):
            position = self.move_down(position)
        return int(position)

    def simulate(self, beads: int, start=None) -> None:
        """
        Show the histogram of results for a specified number of beads on the
        board.

        The simulation condition can be specified by specifying the `start`
        argument (see below).

        Arguments
        ---------
        beads: int
            The number of beads that will be used in the simulation
        start: int, or None
            If `start` is an integer, its value is used by as the starting bin
            position of every bead. If `start` is None (the default), the
            starting bin position of each bead will be randomly chosen from the
            possible bin positions of the board.
        """
        count = Counter(self.run_bead(start) for _ in range(beads))

        (Series(count).sort_index()
                      .mul(100)
                      .div(beads)
                      .plot(kind="bar",
                            xlabel="Bins",
                            ylabel="Relative frequency"))

        plt.xticks(rotation=0)
        plt.gca().yaxis.set_major_formatter(PercentFormatter())
        plt.show()

# Simulate a classic Galton board in which all beads are released at the
# midpoint of the board:
Galton(row_pairs=11, bins=21).simulate(beads=100000, start=11)

# Simulate the variant used in the YouTube show in which any starting position
# is possible:
Galton(row_pairs=11, bins=21).simulate(beads=100000, start=None)

Answer 1

Seems nice.

I'm convinced that the probabilities at the edges are less than the ones in the middle. Mirroring wrong moves using position -= d * 2 will break the uniform distribution. But I guess it may be inside the noise.

Testing random algorithms can be done with dependency injection. You can have some more parameters in __init__ like this:

    def __init__(self, row_pairs: int, bins: int, *, random_position: Callable = random.randint, random_direction: Callable = random.randint):
        (...) 
        self.random_position = random_position
        self.random_direction = random_direction

If you do so, and use these functions in your code, you can inject your own version whenever you want. For example that always give you the least value, or the middle one, or whatever you can implement. With deterministic implementations you can write tests. You can omit the star if you want, but it feels more natural to force these parameters as keyword-only.

Let's say, you want to test what happens if all the beads start from the left, and you always get +1/2. It's one of the easiest cases. You can do it like:

def my_random_position(start, stop):
    return start

def my_random_direction(start, stop):
   return stop

Galton(row_pairs=11, bins=21).simulate(beads=100000, start=11, random_position=my_random_position, random_direction=my_random_direction)

If you also split simulating (creating) data from plotting it. You can write unit tests easily.

You can also use mocking with the same effect, but for me dependency injection is a less magic, more functional, nicer way.

UPDATE:

I would use random.choice([0.5, -0.5]) instead of -0.5 if randint(0, 1) else 0.5

Maybe fractions.Fraction(1/2) would be better however. Using floating point numbers can be problematic if you want to calculate integers because of rounding errors. But 1/2 has a pretty simple binary form, and you not multiply/divide it, so it seems OK.

A started to refactor on Github and I also added some test cases to show what I meant earlier.