2020-03-14 06:51:30 +00:00
|
|
|
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
|
|
|
|
"""
|
2022-01-30 19:29:54 +00:00
|
|
|
return (self.x**2 + self.y**2) <= 1
|
2020-03-14 06:51:30 +00:00
|
|
|
|
|
|
|
@classmethod
|
|
|
|
def random_unit_square(cls):
|
|
|
|
"""
|
|
|
|
Generates a point randomly drawn from the unit square [0, 1) x [0, 1).
|
|
|
|
"""
|
2020-03-14 22:55:13 +00:00
|
|
|
return cls(x=random.random(), y=random.random())
|
|
|
|
|
2020-03-14 06:51:30 +00:00
|
|
|
|
|
|
|
def estimate_pi(number_of_simulations: int) -> float:
|
|
|
|
"""
|
2020-05-01 21:36:35 +00:00
|
|
|
Generates an estimate of the mathematical constant PI.
|
|
|
|
See https://en.wikipedia.org/wiki/Monte_Carlo_method#Overview
|
2020-03-14 06:51:30 +00:00
|
|
|
|
2020-05-01 21:36:35 +00:00
|
|
|
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:
|
2020-03-14 06:51:30 +00:00
|
|
|
|
|
|
|
P[U in unit circle] = 1/4 PI
|
|
|
|
|
|
|
|
and therefore
|
|
|
|
|
|
|
|
PI = 4 * P[U in unit circle]
|
|
|
|
|
2020-05-01 21:36:35 +00:00
|
|
|
We can get an estimate of the probability P[U in unit circle].
|
|
|
|
See https://en.wikipedia.org/wiki/Empirical_probability by:
|
2020-03-14 06:51:30 +00:00
|
|
|
|
|
|
|
1. Draw a point uniformly from the unit square.
|
2020-05-01 21:36:35 +00:00
|
|
|
2. Repeat the first step n times and count the number of points in the unit
|
|
|
|
circle, which is called m.
|
2020-03-14 06:51:30 +00:00
|
|
|
3. An estimate of P[U in unit circle] is m/n
|
|
|
|
"""
|
|
|
|
if number_of_simulations < 1:
|
|
|
|
raise ValueError("At least one simulation is necessary to estimate PI.")
|
|
|
|
|
|
|
|
number_in_unit_circle = 0
|
2022-10-13 16:03:06 +00:00
|
|
|
for _ in range(number_of_simulations):
|
2020-03-14 06:51:30 +00:00
|
|
|
random_point = Point.random_unit_square()
|
|
|
|
|
|
|
|
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()
|
|
|
|
from math import pi
|
2020-03-14 22:55:13 +00:00
|
|
|
|
2020-03-14 06:51:30 +00:00
|
|
|
prompt = "Please enter the desired number of Monte Carlo simulations: "
|
|
|
|
my_pi = estimate_pi(int(input(prompt).strip()))
|
|
|
|
print(f"An estimate of PI is {my_pi} with an error of {abs(my_pi - pi)}")
|