Skip to content

Add doctests and refactor Monte Carlo estimation functions #12703

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Open
wants to merge 2 commits into
base: master
Choose a base branch
from
Open

Add doctests and refactor Monte Carlo estimation functions #12703

Changes from all commits
Commits
Show all changes
2 commits
Select commit Hold shift + click to select a range
bcce276
Refactor Monte Carlo estimation functions and add doctest
akiels Apr 29, 2025
87a5587
cleanup: Remove unused script stub
akiels May 1, 2025
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
192 changes: 117 additions & 75 deletions maths/monte_carlo.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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(
Expand All @@ -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