diff --git a/maths/monte_carlo.py b/maths/monte_carlo.py index d174a0b188a2..af9844c143ec 100644 --- a/maths/monte_carlo.py +++ b/maths/monte_carlo.py @@ -8,35 +8,53 @@ from statistics import mean -def pi_estimator(iterations: int): +def is_in_unit_circle(x: float, y: float) -> bool: """ - An implementation of the Monte Carlo method used to find pi. + Checks whether the point (x, y) lies inside or on the boundary + of a unit circle centered at the origin. + + >>> is_in_unit_circle(0.0, 0.0) + True + >>> is_in_unit_circle(1.0, 0.0) + True + >>> is_in_unit_circle(1.0, 1.0) + False + """ + return sqrt(x**2 + y**2) <= 1 + + +def pi_estimator(iterations: int) -> float: + """ + Estimate the value of pi using the Monte Carlo method. + + The method simulates random points in the square and checks how + many fall within the inscribed unit circle. + + Intuition: 1. Draw a 2x2 square centred at (0,0). - 2. Inscribe a circle within the square. - 3. For each iteration, place a dot anywhere in the square. - a. Record the number of dots within the circle. - 4. After all the dots are placed, divide the dots in the circle by the total. - 5. Multiply this value by 4 to get your estimate of pi. - 6. Print the estimated and numpy value of pi + 2. Inscribe a unit circle within the square. + 3. For each iteration, place a dot randomly within the square. + 4. Divide the number of dots within the circle by the total number of dots. + 5. Multiply this ratio by 4 to estimate pi. + + Args: + iterations: Number of random points to generate. + + Returns: + Estimated value of pi. + + >>> round(pi_estimator(100000), 2) in [3.13, 3.14, 3.15] + True """ - # A local function to see if a dot lands in the circle. - def is_in_circle(x: float, y: float) -> bool: - distance_from_centre = sqrt((x**2) + (y**2)) - # Our circle has a radius of 1, so a distance - # greater than 1 would land outside the circle. - return distance_from_centre <= 1 + if iterations <= 0: + raise ValueError("Number of iterations must be positive") - # The proportion of guesses that landed in the circle - proportion = mean( - int(is_in_circle(uniform(-1.0, 1.0), uniform(-1.0, 1.0))) + inside_circle = sum( + is_in_unit_circle(uniform(-1.0, 1.0), uniform(-1.0, 1.0)) for _ in range(iterations) ) - # The ratio of the area for circle to square is pi/4. - pi_estimate = proportion * 4 - print(f"The estimated value of pi is {pi_estimate}") - print(f"The numpy value of pi is {pi}") - print(f"The total error is {abs(pi - pi_estimate)}") + return (inside_circle / iterations) * 4 def area_under_curve_estimator( @@ -46,86 +64,110 @@ def area_under_curve_estimator( max_value: float = 1.0, ) -> float: """ - An implementation of the Monte Carlo method to find area under - a single variable non-negative real-valued continuous function, - say f(x), where x lies within a continuous bounded interval, - say [min_value, max_value], where min_value and max_value are - finite numbers - 1. Let x be a uniformly distributed random variable between min_value to - max_value - 2. Expected value of f(x) = - (integrate f(x) from min_value to max_value)/(max_value - min_value) - 3. Finding expected value of f(x): - a. Repeatedly draw x from uniform distribution - b. Evaluate f(x) at each of the drawn x values - c. Expected value = average of the function evaluations - 4. Estimated value of integral = Expected value * (max_value - min_value) - 5. Returns estimated value + Estimate the definite integral of a continuous, non-negative function f(x) + over a bounded interval [min_value, max_value] using the Monte Carlo method. + + Intuition: + 1. Randomly sample x values uniformly from [min_value, max_value]. + 2. Compute the average of f(x) at the sample points. + This estimates the expected value of f(x). + 3. Multiply this average by the interval width to estimate the integral. + + Args: + iterations: Number of random samples to draw. + function_to_integrate: A function f(x) to integrate over [min_value, max_value]. + min_value: Lower bound of the interval. + max value: Upper bound of the interval. + + Returns: + Estimated area under the curve (integral of f from min_value to max_value). + + Raises: + ValueError: If iterations <= 0 or min_value >= max_value. + + Example: + >>> def f(x): return x + >>> round(area_under_curve_estimator(100000, f), 2) in [0.49, 0.5, 0.51] + True """ + if iterations <= 0: + raise ValueError("Number of iterations must be positive") + if min_value >= max_value: + raise ValueError("min_value must be less than max_value") - return mean( + samples = ( function_to_integrate(uniform(min_value, max_value)) for _ in range(iterations) - ) * (max_value - min_value) + ) + return mean(samples) * (max_value - min_value) def area_under_line_estimator_check( iterations: int, min_value: float = 0.0, max_value: float = 1.0 -) -> None: +) -> tuple[float, float, float]: """ - Checks estimation error for area_under_curve_estimator function - for f(x) = x where x lies within min_value to max_value - 1. Calls "area_under_curve_estimator" function - 2. Compares with the expected value - 3. Prints estimated, expected and error value + Demonstrate the Monte Carlo integration method by estimating + the area under the line y = x over [min_value, max_value]. + + The exact integral of f(x) = x over [a, b] is: + (b**2 - a**2) / 2 + + This helps illustrate that the Monte Carlo method produces + a reasonable estimate for a simple, well-understood function. + + Args: + iterations: Number of samples to use in estimation. + min_value: Lower bound of the integration interval. + max_value: Upper bound of the integration interval. + + Returns: + A tuple of (estimated_value, expected_value, absolute_error). + + >>> result = area_under_line_estimator_check(100000) + >>> round(result[0], 2) in [0.49, 0.5, 0.51] + True """ def identity_function(x: float) -> float: - """ - Represents identity function - >>> [function_to_integrate(x) for x in [-2.0, -1.0, 0.0, 1.0, 2.0]] - [-2.0, -1.0, 0.0, 1.0, 2.0] - """ return x estimated_value = area_under_curve_estimator( iterations, identity_function, min_value, max_value ) - expected_value = (max_value * max_value - min_value * min_value) / 2 - - print("******************") - print(f"Estimating area under y=x where x varies from {min_value} to {max_value}") - print(f"Estimated value is {estimated_value}") - print(f"Expected value is {expected_value}") - print(f"Total error is {abs(estimated_value - expected_value)}") - print("******************") + expected_value = (max_value**2 - min_value**2) / 2 + error = abs(estimated_value - expected_value) + return estimated_value, expected_value, error -def pi_estimator_using_area_under_curve(iterations: int) -> None: +def pi_estimator_using_area_under_curve(iterations: int) -> tuple[float, float, float]: """ - Area under curve y = sqrt(4 - x^2) where x lies in 0 to 2 is equal to pi + Estimate the value of pi using Monte Carlo integration. + + The area under the curve y = sqrt(4 - x**2) from x = 0 to 2 is equal to pi. + This function uses that fact to estimate pi. + + Args: + iterations: Number of random samples to use for estimation. + + Returns: + A tuple of (estimated_value, expected_value, absolute_error). + + >>> result = pi_estimator_using_area_under_curve(100000) + >>> round(result[0], 2) in [3.13, 3.14, 3.15] + True """ - def function_to_integrate(x: float) -> float: + def semicircle_function(x: float) -> float: """ Represents semi-circle with radius 2 - >>> [function_to_integrate(x) for x in [-2.0, 0.0, 2.0]] + >>> [semicircle_function(x) for x in [-2.0, 0.0, 2.0]] [0.0, 2.0, 0.0] """ return sqrt(4.0 - x * x) estimated_value = area_under_curve_estimator( - iterations, function_to_integrate, 0.0, 2.0 + iterations, semicircle_function, 0.0, 2.0 ) + expected_value = pi + error = abs(estimated_value - expected_value) - print("******************") - print("Estimating pi using area_under_curve_estimator") - print(f"Estimated value is {estimated_value}") - print(f"Expected value is {pi}") - print(f"Total error is {abs(estimated_value - pi)}") - print("******************") - - -if __name__ == "__main__": - import doctest - - doctest.testmod() + return estimated_value, expected_value, error