pi_monte_carlo_estimation.py

import random


class Point:
    def __init__(self, x: float, y: float) -> None:
        self.x = x
        self.y = y

    def is_in_unit_circle(self) -> bool:
        """
        True, if the point lies in the unit circle
        False, otherwise
        """

        return (self.x**2 + self.y**2) <= 1

    @classmethod
    def random_unit_square(cls, ran: random.Random):
        """
        Generates a point randomly drawn from the unit square [0, 1) x [0, 1),
        using 'ran' random number generator
        """

        return cls(x=ran.random(), y=ran.random())


def estimate_pi(number_of_simulations: int, seed: int | None = None) -> float:
    """
    Generates an estimate of the mathematical constant PI.
    See https://en.wikipedia.org/wiki/Monte_Carlo_method#Overview

    The estimate is generated by Monte Carlo simulations. Let U be uniformly drawn from
    the unit square [0, 1) x [0, 1). The probability that U lies in the unit circle is:

        P[U in unit circle] = 1/4 PI

    and therefore

        PI = 4 * P[U in unit circle]

    We can get an estimate of the probability P[U in unit circle].
    See https://en.wikipedia.org/wiki/Empirical_probability by:

        1. Draw a point uniformly from the unit square.
        2. Repeat the first step n times and count the number of points in the unit
            circle, which is called m.
        3. An estimate of P[U in unit circle] is m/n

    'seed' provides seed for the number generator - if None, no seed is used.

    >>> estimate_pi(100, 1)
    3.2
    >>> estimate_pi(1000, 11)
    3.156
    >>> estimate_pi(1000000, 111)
    3.139892

    """
    if number_of_simulations < 1:
        raise ValueError("At least one simulation is necessary to estimate PI.")

    number_in_unit_circle = 0
    ran = random.Random(seed)
    for _ in range(number_of_simulations):
        random_point = Point.random_unit_square(ran)

        if random_point.is_in_unit_circle():
            number_in_unit_circle += 1

    return 4 * number_in_unit_circle / number_of_simulations


if __name__ == "__main__":
    import doctest

    doctest.testmod()