The graphic below shows a square with sides of length 2 containing a circle with a radius of 1. The area of the square is 4 and the area of the circle is 3.141592654 (sound familiar?!). The ratio of the area of the circle to the area of the square is 0.785398163 which is π / 4.
A while ago I wrote an article on estimating pi using various formulae, and in this post I will use the circle-within-a-square ratio illustrated above to run a stochastic simulation in Python to estimate a value of π.
Stochastics
The word stochastic comes from the Greek meaning "to aim at a target", but in general modern usage refers to simulations using random variables.
Such simulations are widely used in many fields from economics to meteorology, including many areas which fall into the category of artificial intelligence and machine learning. The general principle behind stochastics is that it is impossible to accurately capture and model complex systems with vast numbers of data points, but using appropriate random data we can model a simulated system which, by the "law of averages" will closely match the real world. If it does not then the model and/or method of generating random values can be tweaked.
For this project I will combine the "aim at a target" and "simulation using random variables" meanings of stochastic to simulate firing a large number of arrows at the square in the graphic above. If we keep a count of the proportion of arrows which land in the circle this should be somewhere in the region of 0.785398163. When our virtual quiver is empty we can multiply the ratio we end up with by 4 to get an estimation of π.
This is not of course the best way of estimating pi and using it to do so is frankly overkill. However, it does give a simple introduction into the principle of using random variables in a simulation.
Coding
Create a folder somewhere handy and within it create the following empty files. You can download them as a zip or clone/download the Github repository if you prefer.
- stochasticpi.py
- main.py
Source Code Links
Open stochasticpi.py and type or paste the following. You can download the files in a zip from the Downloads page or clone/download from Github if you prefer.
stochasticpi.py
import random import math def start(arrow_count, output_interval): """ Simulate firing arrows at a circle within a square, the proportion hitting the circle being used to estimate a value for Pi. """ arrows_fired = 0 circle_hits = 0 pi_estimate = 0 x = 0 y = 0 for i in range(0, arrow_count): x = random.uniform(-1, 1) y = random.uniform(-1, 1) arrows_fired += 1 if _in_unit_circle(x, y): circle_hits += 1 pi_estimate = (circle_hits / arrows_fired) * 4 if arrows_fired % output_interval == 0: print("Intermediate estimate of Pi after {} arrows: {:.16f}".format(arrows_fired, pi_estimate)) print("Final estimate of Pi after {} arrows: {:.16f}".format(arrow_count, pi_estimate)) def _in_unit_circle(x, y): """ Calculate whether a particular position within a square is within a circle fitting the square. """ if math.sqrt(x**2 + y**2) <= 1: return True else: return False
At the top of the file we import random and math, and then within start we first set up a few variables, all initialized to 0. We then enter a loop from 0 to the number of iterations. Firing an arrow is simulated by setting x and y to random values between -1 and 1. (This assumes the centre of the square is 0,0.) We then call a separate function to check if the arrow hit the circle; if so circle_hits is incremented. We then update the estimate of pi. At the end of the loop we check if the current iteration is a multiple of the output interval; if so we print the estimate arrived at so far. This process slows everything down a lot so if you are just interested in the end result you could remove it as well as the output_interval argument. If you do then you will only need to calculate pi_estimate once at the end.
The process of checking whether an arrow has hit the circle is carried out by a separate function called _in_unit_circle. This uses the old "square on the hypotenuse is equal to...", the values of x and y being "the other two sides" and the line from the centre of the circle to x,y being the hypotenuse. If the hypotenuse is <= 1 then the point is inside the circle, otherwise it is outside.
That's stochasticpi.py finished so let's write a main function to try it out.
main.py
import stochasticpi def main(): """ Here we start the stochastic Pi simulation. Arguments are total number of iterations and interval between showing results. """ print("-----------------") print("| codedrome.com |") print("| Stochastic Pi |") print("-----------------\n") stochasticpi.start(10000000, 1000000) main()
Just a single line in main to call the start function, so run the program with this command:
Running the program
python3.7 main.py
The output is:
Program Output
----------------- | codedrome.com | | Stochastic Pi | ----------------- Intermediate estimate of Pi after 1000000 arrows: 3.1415560000000000 Intermediate estimate of Pi after 2000000 arrows: 3.1416279999999999 Intermediate estimate of Pi after 3000000 arrows: 3.1424853333333331 Intermediate estimate of Pi after 4000000 arrows: 3.1419850000000000 Intermediate estimate of Pi after 5000000 arrows: 3.1416784000000000 Intermediate estimate of Pi after 6000000 arrows: 3.1421006666666669 Intermediate estimate of Pi after 7000000 arrows: 3.1419268571428574 Intermediate estimate of Pi after 8000000 arrows: 3.1418119999999998 Intermediate estimate of Pi after 9000000 arrows: 3.1418320000000000 Intermediate estimate of Pi after 10000000 arrows: 3.1417028000000000 Final estimate of Pi after 10000000 arrows: 3.1417028000000000
Every run will be different and in this run I only got 3dp of accuracy despite using a million arrows. However, the purpose of this post was to demonstrate a technique rather than develop a practical solution.