From cdbd8a5e3871b734c362153342251450dab504f1 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Christian=20Sch=C3=BCrch?= Date: Thu, 23 Jan 2020 16:25:46 +0100 Subject: [PATCH 1/5] Add Monte Carlo estimation of PI --- maths/pi_monte_carlo_estimation.py | 63 ++++++++++++++++++++++++++++++ 1 file changed, 63 insertions(+) create mode 100644 maths/pi_monte_carlo_estimation.py diff --git a/maths/pi_monte_carlo_estimation.py b/maths/pi_monte_carlo_estimation.py new file mode 100644 index 000000000000..d6367f81da1a --- /dev/null +++ b/maths/pi_monte_carlo_estimation.py @@ -0,0 +1,63 @@ +import random + + +class Point: + def __init__(self, x, y): + self.x = x + self.y = y + + def is_in_unit_circle(self): + """ + 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): + """ + Generates a point randomly drawn from the unit square [0, 1) x [0, 1). + """ + x = random.random() + y = random.random() + + return cls(x, y) + + +def estimate_pi(number_of_simulations): + """ + 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 + """ + if number_of_simulations < 1: + raise ValueError("At least one simulation is necessary to estimate PI.") + + number_in_unit_circle = 0 + for simulation_index in range(number_of_simulations): + 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__": + number_of_simulations = int( + input("Please enter the number of Monte Carlo simulations: ").strip() + ) + + print(f"An estimate of PI is: {estimate_pi(number_of_simulations)}") From eaff1422d5d7eb0e58695a85f89fa472e13cb4ac Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Christian=20Sch=C3=BCrch?= Date: Fri, 24 Jan 2020 11:28:47 +0100 Subject: [PATCH 2/5] Add type annotations for Monte Carlo estimation of PI --- maths/pi_monte_carlo_estimation.py | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/maths/pi_monte_carlo_estimation.py b/maths/pi_monte_carlo_estimation.py index d6367f81da1a..b4d3d55fd25d 100644 --- a/maths/pi_monte_carlo_estimation.py +++ b/maths/pi_monte_carlo_estimation.py @@ -2,16 +2,16 @@ class Point: - def __init__(self, x, y): + def __init__(self, x: float, y: float) -> None: self.x = x self.y = y - def is_in_unit_circle(self): + 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 + return (self.x ** 2 + self.y ** 2) <= 1 @classmethod def random_unit_square(cls): @@ -24,7 +24,7 @@ def random_unit_square(cls): return cls(x, y) -def estimate_pi(number_of_simulations): +def estimate_pi(number_of_simulations: int) -> float: """ Generates an estimate of the mathematical constant PI (see https://en.wikipedia.org/wiki/Monte_Carlo_method#Overview). From 67bcc8dde9e8f5d2361b6578c12866193592e46a Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Christian=20Sch=C3=BCrch?= Date: Fri, 24 Jan 2020 12:45:45 +0100 Subject: [PATCH 3/5] Compare the PI estimate to PI from the math lib --- maths/pi_monte_carlo_estimation.py | 10 ++++++---- 1 file changed, 6 insertions(+), 4 deletions(-) diff --git a/maths/pi_monte_carlo_estimation.py b/maths/pi_monte_carlo_estimation.py index b4d3d55fd25d..fca5b27b6b17 100644 --- a/maths/pi_monte_carlo_estimation.py +++ b/maths/pi_monte_carlo_estimation.py @@ -56,8 +56,10 @@ def estimate_pi(number_of_simulations: int) -> float: if __name__ == "__main__": - number_of_simulations = int( - input("Please enter the number of Monte Carlo simulations: ").strip() - ) + # import doctest - print(f"An estimate of PI is: {estimate_pi(number_of_simulations)}") + # doctest.testmod() + from math import pi + 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 accuracy of {abs(my_pi - pi)}") From 0a990ac597d9bde06326a49372cd52d70a170d49 Mon Sep 17 00:00:00 2001 From: John Law Date: Sat, 14 Mar 2020 06:49:17 +0100 Subject: [PATCH 4/5] accuracy -> error --- maths/pi_monte_carlo_estimation.py | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/maths/pi_monte_carlo_estimation.py b/maths/pi_monte_carlo_estimation.py index fca5b27b6b17..43d4ac7cec84 100644 --- a/maths/pi_monte_carlo_estimation.py +++ b/maths/pi_monte_carlo_estimation.py @@ -23,7 +23,6 @@ def random_unit_square(cls): return cls(x, y) - def estimate_pi(number_of_simulations: int) -> float: """ Generates an estimate of the mathematical constant PI (see https://en.wikipedia.org/wiki/Monte_Carlo_method#Overview). @@ -62,4 +61,4 @@ def estimate_pi(number_of_simulations: int) -> float: from math import pi 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 accuracy of {abs(my_pi - pi)}") + print(f"An estimate of PI is {my_pi} with an error of {abs(my_pi - pi)}") From 18932bcdeb53f42511d8251d99fb64794e859d88 Mon Sep 17 00:00:00 2001 From: John Law Date: Sat, 14 Mar 2020 07:44:37 +0100 Subject: [PATCH 5/5] Update pi_monte_carlo_estimation.py --- maths/pi_monte_carlo_estimation.py | 5 +---- 1 file changed, 1 insertion(+), 4 deletions(-) diff --git a/maths/pi_monte_carlo_estimation.py b/maths/pi_monte_carlo_estimation.py index 43d4ac7cec84..7f341ade94a4 100644 --- a/maths/pi_monte_carlo_estimation.py +++ b/maths/pi_monte_carlo_estimation.py @@ -18,10 +18,7 @@ def random_unit_square(cls): """ Generates a point randomly drawn from the unit square [0, 1) x [0, 1). """ - x = random.random() - y = random.random() - - return cls(x, y) + return cls(x = random.random(), y = random.random()) def estimate_pi(number_of_simulations: int) -> float: """