From 1bc84e1fa0b3bbd75c73b0b1b25f573c241e1c9b Mon Sep 17 00:00:00 2001
From: cschuerc <57899042+cschuerc@users.noreply.github.com>
Date: Sat, 14 Mar 2020 07:51:30 +0100
Subject: [PATCH] Add Monte Carlo estimation of PI (#1712)

* Add Monte Carlo estimation of PI

* Add type annotations for Monte Carlo estimation of PI

* Compare the PI estimate to PI from the math lib

* accuracy -> error

* Update pi_monte_carlo_estimation.py

Co-authored-by: John Law <johnlaw.po@gmail.com>
---
 maths/pi_monte_carlo_estimation.py | 61 ++++++++++++++++++++++++++++++
 1 file changed, 61 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 000000000..7f341ade9
--- /dev/null
+++ b/maths/pi_monte_carlo_estimation.py
@@ -0,0 +1,61 @@
+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):
+        """
+        Generates a point randomly drawn from the unit square [0, 1) x [0, 1).
+        """
+        return cls(x = random.random(), y = random.random())
+
+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).
+
+    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__":
+    # import doctest
+
+    # 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 error of {abs(my_pi - pi)}")