Description
1. The idea of using a probabilistic (random) simulation to model a physical process is called a Monte Carlo simulation.
We can simulate a physical process to create an estimate of π. The physical process in question is that of throwing darts at a target. Imagine a square board of side 2 units, with a circular target of radius 1 unit drawn in the centre (0,0) of the board:
The idea: Throw a large (known) number of darts at the board. The fraction of the darts that fall within the circle will be proportional to the area of the circle. Assume all darts thrown hit the board:
Area of Circle π(1)2 Num. darts in circle
= =
Area of Square (2)2 Total num. darts
The mechanics of doing the simulation simplify considerably if we conduct the simulation in the first quadrant of the target:
The part of the square has corners (0, 0) (1,0) (1,1) (0,1). The quarter circle has radius 1 unit. If we throw darts at the quarter square then, following the calculation above:
Area of 1/4 Circle π(1)2/4 Num. darts in 1/4 circle
= =
Area of 1/4 Square (1)2 Total num. darts thrown
From this, we obtain an estimate of π:
4 × Num darts in 1/4 circle πest = . Total num. Darts thrown
• Throwing a dart: A thrown dart falls (randomly) on the 1/4 board.
The coordinates of the dart hits are: 0 ≤ xthrow < 1 and 0 ≤ ythrow < 1. Each throw can be simulated by a random number generator. Obtain this by importing the random module. Then x throw = random.random() generates a random float in between 0 and 1 and similarly for y throw
• Count hits: To count darts that hit the circle, we use the throw above, with coordinates (xthrow,ythrow) and check that it falls within a circle of radius r = 1.
(hint: the boundary of the red circle is given by x2 + y2 = r2 = (1)2.) All dart throws: (xi,yi) (1 ≤ i ≤ N) such that x2i + yi2 ≤ r2 = 1, fall inside the circle.
• Run simulations of N = 10,100,1000,10000 dart throws and compare the estimates of π. Use math.pi as a nominative “exact” value. Make a table: (Number of trials, Estimate, Relative Error). Note: Relative Error = (measured – exact)/exact
Storing results in a numpy array.
A one-dimensional numpy array is a list of values, all of the same type, indexed by nonnegative integers.
• Create an numpy array of size about n = 6.
• (import numpy as np at the top of the program).
• np.zeros(k) produces [0., 0., …, 0] (an array of k zeros)
(It also serves to initialize the array.)
• Store the estimates, by range, in the array.
For example, a[0] = 7 assigns 7 to the first array element.
Suppose most of the estimates obtained in the simulations are in the range: (3.10, 3.20), with a few outside that range.
• Store the count of estimates that are less than 3.10 in array loc 0.
• Store estimates greater than 3.10 & less than 3.12 in array loc.1 etc.
• Each simulation should have 10,000 dart throws.
• Run 10,000 such simulations. Print out the number of estimates of π in each range.
2. Making a histogram
This exercise makes a histogram from the data gathered in the simulation above. There you ran 10,000 simulations of 10,000 dart throws each, and stored the estimates in an array. Use the data generated by your program to construct a histogram of the 10,000 estimates in your array: Use Matplotlib for making the histogram. To use Matplotlib:
First import matplotlib.pyplot as plt, the command
plt.hist(x, nBins) makes the histogram, where x is the array containing the data and nBins is a variable setting the number of bins.
Use: plt.xlabel(’ π Estimates ’) plt.ylabel(’ Counts ’) plt.title(’Histogram of π estimates’)
https://matplotlib.org/api/_as_gen/matplotlib.pyplot.hist.html The above documentation contains examples in a gallery at the bottom of the page.
An alternate way to make a histogtram:
(i) Create the individual bars using: fig, bars = plt.subplots()
(ii) Let the count data (i.e. the number of values in each interval) be in an array y. Create the bins using:
x = np.arange(len(y)) ( np.arange(stop) acts just like range(stop).) https://docs.scipy.org/doc/numpy/reference/generated/numpy.arange.html# numpy.arange)
For either method, plt.show() should be the final command.
3. Write a script to simulate the following experiment: Mollie and Chloe agree to meet at the coffee shop between 9:30 and 10:00 am. They are pretty casual about the logistics. Each one will show up at some time during the half hour interval. If Mollie arrives first, she will wait 5 minutes and leave if Chloe does not arrive by then. Similarly, if Chloe arrives first, she will wait 7 minutes for Mollie. Neither waits past 10:00am.
• What is the probability that they meet?
• How does the answer change if Chloe reduces her waiting time to 5 to match Mollie’s?
• How does the answer change if Mollie increases her waiting time to 7 to match Chloe’s?
