EGR 103/Fall 2018/Lec 12

This page will go through the development of the program we wrote in class. Each step along the way, I will describe what functionality we wanted to add and you will be able to see the code that accomplished that goal.

Introduction

We will be writing a program that uses a Monte Carlo method to estimate \(\pi\). The program will simulate throwing some number of darts at a dart board of radius \(r\) mounted on a square back board that is \(2r\) on a side. If the \(x\) and </math>y</math> coordinates are each uniformly distributed between \(-r\) and \(r\), then the probability that a dart hits the dart board is the same as the ratio of the area of the dart board to the area of the back board - that is, \(p_{hit}=\frac{\pi r^2}{(2r)^2}=\frac{\pi}{4}\). We will thus estimate \(\pi\) by dividing the number of darts that hit the dart board by the total number of darts thrown and then multiplying that by 4.

Detached figure windows

As a part of this code, we will be watching a plot as it evolves. That is not something you want to do with inline graphics. To make Spyder put figures in their own windows:

• For Windows, go to Tools->Preferences and with OSX, go to python->Preferences
• In the Preferences window, select IPython console from the list at left
• In the right window, select the Graphics tab
• In the Graphics tab, in the Graphics backend section, choose Automatic.
• Click apply, then get out of the Preferences window. You will need to restart Spyder for this to take effect.

Tasks

We broke the tasks up into several parts and wrote to code to accomplish those tasks. The main tasks were:

• Import modules
• Draw the board
• Throw a dart
• Throw lots of darts
• Plot the dart locations
• Determine how many darts land on the dart board
• Print how many darts land on the dart board and the estimates of \(\pi\)
• Refine how that is printed
• Plot the evolution of the estimates
• Have the value of a variable decide if we should plot while the simulation is running or only at the end
• Figure out what's next

Importing modules

As always, start by importing the necessary modules. We will be doing math, using computational methods, and making graphs, so we will need:

1 import math as m
2 import numpy as np
3 import matplotlib.pyplot as plt

Drawing the board

Assuming a dart board of radius \(r\), the square that encloses it would be \(2r\) on a side. Since we will later be using the distance from the center of the circle to determine if a dart is on the dart board, it is mathematically easiest to have the middle of the board at (0, 0). Drawing the back board thus means drawing a square with vertices at (-r, -r), (r, -r), (r, r), and (-r, r). The the circle, we will need a list of x and y coordinates that defines a circle; the easiest way to do that is to use polar coordinates with a constant radius and an angle that goes from 0 to \(2\pi\). We need to know how many values to use in making the circle, so that will be another input to the function. We chose a default radius of 1 and a default number of points of 100:

 1 import math as m
 2 import numpy as np
 3 import matplotlib.pyplot as plt
 4 
 5 def draw_board(r=1, pts=100):
 6     plt.plot([-r, r, r, -r, -r], [-r, -r, r, r, -r], 'k-')
 7     theta = np.linspace(0, 2*m.pi, pts)
 8     x = r*np.cos(theta)
 9     y = r*np.sin(theta)
10     plt.plot(x, y, 'r-')

Test the code by calling it from the script:

 1 import math as m
 2 import numpy as np
 3 import matplotlib.pyplot as plt
 4 
 5 def draw_board(r=1, pts=100):
 6     plt.plot([-r, r, r, -r, -r], [-r, -r, r, r, -r], 'k-')
 7     theta = np.linspace(0, 2*m.pi, pts)
 8     x = r*np.cos(theta)
 9     y = r*np.sin(theta)
10     plt.plot(x, y, 'r-')
11 
12 if __name__ == '__main__':
13     plt.figure(1)
14     plt.clf()
15     draw_board()

The problem with this is that Python will use all the figure window and depending on the shape of the window, the square and the circle might look non-square and non-circular. To require that Python plot using the same scale along both axes, use the .axis('equal') method of plt:

 1 import math as m
 2 import numpy as np
 3 import matplotlib.pyplot as plt
 4 
 5 def draw_board(r=1, pts=100):
 6     plt.plot([-r, r, r, -r, -r], [-r, -r, r, r, -r], 'k-')
 7     theta = np.linspace(0, 2*m.pi, pts)
 8     x = r*np.cos(theta)
 9     y = r*np.sin(theta)
10     plt.plot(x, y, 'r-')
11     plt.axis('equal')
12 
13 if __name__ == '__main__':
14     plt.figure(1)
15     plt.clf()
16     draw_board()

You will now have a square back board and a circular dart board. As an aside, if you want to see something trippy, copy the following into the console after you have run your script:

r = 1
for k in range(10):
     draw_board(r)
     r = r*m.sqrt(2)

Throw a dart

Throwing a dart means figuring out the \(x\) and \(y\) coordinates for a dart. To determine the valid location of darts, we need to know the radius of the dart board. Given all that, adding a function to return the location of a thrown dart might look like:

 1 import math as m
 2 import numpy as np
 3 import matplotlib.pyplot as plt
 4 
 5 def draw_board(r=1, pts=100):
 6     plt.plot([-r, r, r, -r, -r], [-r, -r, r, r, -r], 'k-')
 7     theta = np.linspace(0, 2*m.pi, pts)
 8     x = r*np.cos(theta)
 9     y = r*np.sin(theta)
10     plt.plot(x, y, 'r-')
11     plt.axis('equal')
12     
13 def throw_darts(r=1):
14     x = np.random.uniform(-r, r)
15     y = np.random.uniform(-r, r)
16     return (x, y)
17 
18 if __name__ == '__main__':
19     plt.figure(1)
20     plt.clf()
21     draw_board()

You can test this in the console by throwing darts at a board of a given size (say, 10):

for k in range(10):
    print(throw_darts(10))

Make sure both coordinates are in the proper range.

Throw lots of darts

Now that you can throw a single dart, you might want to throw several darts. You should add variables to keep track of the total number of darts to throw and of the total number thrown:

 1 import math as m
 2 import numpy as np
 3 import matplotlib.pyplot as plt
 4 
 5 def draw_board(r=1, pts=100):
 6     plt.plot([-r, r, r, -r, -r], [-r, -r, r, r, -r], 'k-')
 7     theta = np.linspace(0, 2*m.pi, pts)
 8     x = r*np.cos(theta)
 9     y = r*np.sin(theta)
10     plt.plot(x, y, 'r-')
11     plt.axis('equal')
12     
13 def throw_darts(r=1):
14     x = np.random.uniform(-r, r)
15     y = np.random.uniform(-r, r)
16     return (x, y)
17     
18 if __name__ == '__main__':
19     darts_to_throw = 100
20     darts_thrown = 0
21     
22     plt.figure(1)
23     plt.clf()
24     draw_board()
25 
26     while darts_thrown < darts_to_throw:
27         darts_thrown += 1
28         x_d, y_d = throw_darts()

Plot the dart locations

Now that you are throwing darts, you can plot the locations at which they land. There are two fundamentally different ways to do this - you can plot each dart as it is thrown or you can wait until the end to plot all the darts. We are going to first look at how to plot each dart and watch the board fill up with darts. They key here will be a function that takes the dart coordinates and then updates the figure before moving on to the next command. Having the figure update right now requires drawing the plot and then adding a pause to let the plot...plot:

 1 import math as m
 2 import numpy as np
 3 import matplotlib.pyplot as plt
 4 
 5 def draw_board(r=1, pts=100):
 6     plt.plot([-r, r, r, -r, -r], [-r, -r, r, r, -r], 'k-')
 7     theta = np.linspace(0, 2*m.pi, pts)
 8     x = r*np.cos(theta)
 9     y = r*np.sin(theta)
10     plt.plot(x, y, 'r-')
11     plt.axis('equal')
12     
13 def throw_darts(r=1):
14     x = np.random.uniform(-r, r)
15     y = np.random.uniform(-r, r)
16     return (x, y)
17     
18 def plot_darts(x, y):
19     plt.figure(1)
20     plt.plot(x, y, 'b.')
21     plt.draw()
22     plt.pause(0.01)
23     
24 if __name__ == '__main__':
25     darts_to_throw = 100
26     darts_thrown = 0
27     
28     plt.figure(1)
29     plt.clf()
30     draw_board()
31 
32     while darts_thrown < darts_to_throw:
33         darts_thrown += 1
34         x_d, y_d = throw_darts()
35         plot_darts(x_d, y_d)

Determine how many darts land on the dart board

Now that you are throwing darts and see where they hit, you can add code to keep track of how many are landing on the dart board. To do this, you will need to let your code know how big the dart board is and then calculate how far the throw is from the center. You need to keep track of how many darts hit the dartboard.

 1 import math as m
 2 import numpy as np
 3 import matplotlib.pyplot as plt
 4 
 5 def draw_board(r=1, pts=100):
 6     plt.plot([-r, r, r, -r, -r], [-r, -r, r, r, -r], 'k-')
 7     theta = np.linspace(0, 2*m.pi, pts)
 8     x = r*np.cos(theta)
 9     y = r*np.sin(theta)
10     plt.plot(x, y, 'r-')
11     plt.axis('equal')
12     
13 def throw_darts(r=1):
14     x = np.random.uniform(-r, r)
15     y = np.random.uniform(-r, r)
16     return (x, y)
17     
18 def plot_darts(x, y):
19     plt.figure(1)
20     plt.plot(x, y, 'b.')
21     plt.draw()
22     plt.pause(0.01)
23     
24 if __name__ == '__main__':
25     darts_to_throw = 100
26     darts_thrown = 0
27     board_rad = 1
28     inside = 0
29 
30     plt.figure(1)
31     plt.clf()
32     draw_board(r=board_rad)
33 
34     while darts_thrown < darts_to_throw:
35         darts_thrown += 1
36         x_d, y_d = throw_darts(board_rad)
37         plot_darts(x_d, y_d)
38         r_d = m.sqrt(x_d**2 + y_d**2)
39         if r_d < board_rad:
40             inside += 1

Note the changes on lines 32 and 36 to allow for board sizes other than 1.

Print how many darts land on the dart board and the estimates of \(\pi\)

While you can see the value of inside in the variable explorer, it would make more sense to print information out as the simulation is running. This might include how many darts have been thrown, how many hit the dart board, and what the estimate for \(\pi\) is:

 1 import math as m
 2 import numpy as np
 3 import matplotlib.pyplot as plt
 4 
 5 def draw_board(r=1, pts=100):
 6     plt.plot([-r, r, r, -r, -r], [-r, -r, r, r, -r], 'k-')
 7     theta = np.linspace(0, 2*m.pi, pts)
 8     x = r*np.cos(theta)
 9     y = r*np.sin(theta)
10     plt.plot(x, y, 'r-')
11     plt.axis('equal')
12     
13 def throw_darts(r=1):
14     x = np.random.uniform(-r, r)
15     y = np.random.uniform(-r, r)
16     return (x, y)
17     
18 def plot_darts(x, y):
19     plt.figure(1)
20     plt.plot(x, y, 'b.')
21     plt.draw()
22     plt.pause(0.01)
23     
24 if __name__ == '__main__':
25     darts_to_throw = 100
26     darts_thrown = 0
27     board_rad = 1
28     inside = 0
29 
30     plt.figure(1)
31     plt.clf()
32     draw_board(r=board_rad)
33 
34     while darts_thrown < darts_to_throw:
35         darts_thrown += 1
36         x_d, y_d = throw_darts(board_rad)
37         plot_darts(x_d, y_d)
38         r_d = m.sqrt(x_d**2 + y_d**2)
39         if r_d < board_rad:
40             inside += 1
41 
42         print(darts_thrown, inside, 4*inside/darts_thrown)

Refine how that is printed

That table is fine for a rough idea of how things are going, but it would be better to make a formatted table. Assuming you throw fewer than 100,000 darts, you can format the number of darts thrown and the number of hits to be integers with 5 spaces reserved and the estimate of pi as a floating point number using scientific notation with three digits after the decimal point. Note the spaces in the format strings - those make sure the three numbers are separated.

 1 import math as m
 2 import numpy as np
 3 import matplotlib.pyplot as plt
 4 
 5 def draw_board(r=1, pts=100):
 6     plt.plot([-r, r, r, -r, -r], [-r, -r, r, r, -r], 'k-')
 7     theta = np.linspace(0, 2*m.pi, pts)
 8     x = r*np.cos(theta)
 9     y = r*np.sin(theta)
10     plt.plot(x, y, 'r-')
11     plt.axis('equal')
12     
13 def throw_darts(r=1):
14     x = np.random.uniform(-r, r)
15     y = np.random.uniform(-r, r)
16     return (x, y)
17     
18 def plot_darts(x, y):
19     plt.figure(1)
20     plt.plot(x, y, 'b.')
21     plt.draw()
22     plt.pause(0.01)
23     
24 if __name__ == '__main__':
25     darts_to_throw = 100
26     darts_thrown = 0
27     board_rad = 1
28     inside = 0
29 
30     plt.figure(1)
31     plt.clf()
32     draw_board(r=board_rad)
33 
34     while darts_thrown < darts_to_throw:
35         darts_thrown += 1
36         x_d, y_d = throw_darts(board_rad)
37         plot_darts(x_d, y_d)
38         r_d = m.sqrt(x_d**2 + y_d**2)
39         if r_d < board_rad:
40             inside += 1
41 
42         print('{:5d}'.format(darts_thrown) + 
43                ' {:5d}'.format(inside) +
44                ' {:0.3e}'.format(4*inside/darts_thrown))

Plot the evolution of the estimates

It may be interesting to see how the estimate for \(\pi\) changes as the simulation runs. While the table now prints this information, making a graph might be useful. We can once again either make a graph by adding each estimate or by graphing all the estimates at the end. We will do the former here by creating a second figure window and plotting the estimates on it:

 1 import math as m
 2 import numpy as np
 3 import matplotlib.pyplot as plt
 4 
 5 def draw_board(r=1, pts=100):
 6     plt.plot([-r, r, r, -r, -r], [-r, -r, r, r, -r], 'k-')
 7     theta = np.linspace(0, 2*m.pi, pts)
 8     x = r*np.cos(theta)
 9     y = r*np.sin(theta)
10     plt.plot(x, y, 'r-')
11     plt.axis('equal')
12     
13 def throw_darts(r=1):
14     x = np.random.uniform(-r, r)
15     y = np.random.uniform(-r, r)
16     return (x, y)
17     
18 def plot_darts(x, y):
19     plt.figure(1)
20     plt.plot(x, y, 'b.')
21     plt.draw()
22     plt.pause(0.01)
23 
24 def plot_est(nums, ests):
25     plt.figure(2)
26     plt.plot(nums, ests, 'ko')
27     plt.draw()
28     plt.pause(0.01)
29     
30 if __name__ == '__main__':
31     darts_to_throw = 100
32     darts_thrown = 0
33     board_rad = 1
34     inside = 0
35 
36     plt.figure(1)
37     plt.clf()
38     draw_board(r=board_rad)
39 
40     plt.figure(2)
41     plt.clf()
42 
43     while darts_thrown < darts_to_throw:
44         darts_thrown += 1
45         x_d, y_d = throw_darts(board_rad)
46         plot_darts(x_d, y_d)
47         r_d = m.sqrt(x_d**2 + y_d**2)
48         if r_d < board_rad:
49             inside += 1
50 
51         print('{:5d}'.format(darts_thrown) + 
52                ' {:5d}'.format(inside) +
53                ' {:0.3e}'.format(4*inside/darts_thrown))
54         plot_est(darts_thrown, 4*inside/darts_thrown)

While this works, we are having to calculate the estimates twice. It might make more sense to store all the estimates and then give the most recently calculated estimate to both the line that prints the estimate and the one that plots it:

 1 import math as m
 2 import numpy as np
 3 import matplotlib.pyplot as plt
 4 
 5 def draw_board(r=1, pts=100):
 6     plt.plot([-r, r, r, -r, -r], [-r, -r, r, r, -r], 'k-')
 7     theta = np.linspace(0, 2*m.pi, pts)
 8     x = r*np.cos(theta)
 9     y = r*np.sin(theta)
10     plt.plot(x, y, 'r-')
11     plt.axis('equal')
12     
13 def throw_darts(r=1):
14     x = np.random.uniform(-r, r)
15     y = np.random.uniform(-r, r)
16     return (x, y)
17     
18 def plot_darts(x, y):
19     plt.figure(1)
20     plt.plot(x, y, 'b.')
21     plt.draw()
22     plt.pause(0.01)
23 
24 def plot_est(nums, ests):
25     plt.figure(2)
26     plt.plot(nums, ests, 'ko')
27     plt.draw()
28     plt.pause(0.01)
29     
30 if __name__ == '__main__':
31     darts_to_throw = 100
32     darts_thrown = 0
33     board_rad = 1
34     inside = 0
35     pi_est = []
36 
37     plt.figure(1)
38     plt.clf()
39     draw_board(r=board_rad)
40 
41     plt.figure(2)
42     plt.clf()
43 
44     while darts_thrown < darts_to_throw:
45         darts_thrown += 1
46         x_d, y_d = throw_darts(board_rad)
47         plot_darts(x_d, y_d)
48         r_d = m.sqrt(x_d**2 + y_d**2)
49         if r_d < board_rad:
50             inside += 1
51 
52         pi_est += [4*inside/darts_thrown]
53         print('{:5d}'.format(darts_thrown) + 
54                ' {:5d}'.format(inside) +
55                ' {:0.3e}'.format(pi_est[-1]))
56         plot_est(darts_thrown, pi_est[-1])

Have the value of a variable decide if we should plot while the simulation is running or only at the end

Watching the animation is great, but very slow. If I want to calculate the estimate after 50,000 trials, that will take a great deal of time. Instead, I can add a variable that will indicate to the program whether it should plot as the simulation is running or only at the end. For the latter, we will need to keep track of all the dart throws - not just the current one. The do_plot variable will be the determining factor of whether I plot in the while loop or only when it is done. Pay careful attention to the indentations and also the fact that the plot commands have been moved to the end of the loop:

 1 import math as m
 2 import numpy as np
 3 import matplotlib.pyplot as plt
 4 
 5 def draw_board(r=1, pts=100):
 6     plt.plot([-r, r, r, -r, -r], [-r, -r, r, r, -r], 'k-')
 7     theta = np.linspace(0, 2*m.pi, pts)
 8     x = r*np.cos(theta)
 9     y = r*np.sin(theta)
10     plt.plot(x, y, 'r-')
11     plt.axis('equal')
12     
13 def throw_darts(r=1):
14     x = np.random.uniform(-r, r)
15     y = np.random.uniform(-r, r)
16     return (x, y)
17     
18 def plot_darts(x, y):
19     plt.figure(1)
20     plt.plot(x, y, 'b.')
21     plt.draw()
22     plt.pause(0.01)
23 
24 def plot_est(nums, ests):
25     plt.figure(2)
26     plt.plot(nums, ests, 'ko')
27     plt.draw()
28     plt.pause(0.01)
29     
30 if __name__ == '__main__':
31     darts_to_throw = 100
32     darts_thrown = 0
33     board_rad = 1
34     do_plot = 1
35     inside = 0
36     pi_est = []
37     x_vec = []
38     y_vec = []
39 
40     plt.figure(1)
41     plt.clf()
42     draw_board(r=board_rad)
43 
44     plt.figure(2)
45     plt.clf()
46 
47     while darts_thrown < darts_to_throw:
48         darts_thrown += 1
49         x_d, y_d = throw_darts(board_rad)
50         x_vec += [x_d]
51         y_vec += [y_d]
52         r_d = m.sqrt(x_d**2 + y_d**2)
53         if r_d < board_rad:
54             inside += 1
55 
56         pi_est += [4*inside/darts_thrown]
57         print('{:5d}'.format(darts_thrown) + 
58                ' {:5d}'.format(inside) +
59                ' {:0.3e}'.format(pi_est[-1]))
60         
61         if do_plot:
62             plot_darts(x_d, y_d)
63             plot_est(darts_thrown, pi_est[-1])
64 
65     if not do_plot:
66         plot_darts(x_vec, y_vec)
67         plot_est(np.linspace(1, darts_thrown, darts_thrown), pi_est)

Figure out what's next

The program will nor run through a pre-determined number of simulations, printing a table of estimates as it goes and plotting the evolution of the dart board and the estimates either during or after the simulation. The question is - what's next? There were several recommendations during class, including:

• Re-print the dart board at the end so you can see it even when there are thousands of blue dots - this requires adding:

    plt.figure(1)
    draw_board(board_rad)

to the end of the code

• Have the actual value of \(\pi\) show up on the estimate graph - this requires changing the plot_est function to:

def plot_est(nums, ests):
    plt.figure(2)
    plt.plot(nums, ests, 'ko')
    plt.plot(nums, m.pi*np.ones(np.size(nums)), 'g.')
    plt.draw()
    plt.pause(0.01)

Note that the np.ones(np.size(nums)) is there to make sure there are as many y values as there are x values in case nums is a list.

• Show the darts in different colors - several ways to do this but you have to carefully consider that the plot_darts command is currently set to work with both individual points and lists of points. It would be more complicated to do this with the whole list of points since there is not currently a mechanism to decide which are in and which are out. Some options include:
  • Only choosing different colors when plotting one point at a time and plotting the whole list in a single color
  • Keeping track of which points are inside and which are outside and plotting each list separately - here is an example of the code with changes needed to make that happen (and also re-print the board and show the actual value of pi on the estimate evolution):

 1 import math as m
 2 import numpy as np
 3 import matplotlib.pyplot as plt
 4 
 5 def draw_board(r=1, pts=100):
 6     plt.plot([-r, r, r, -r, -r], [-r, -r, r, r, -r], 'k-')
 7     theta = np.linspace(0, 2*m.pi, pts)
 8     x = r*np.cos(theta)
 9     y = r*np.sin(theta)
10     plt.plot(x, y, 'r-')
11     plt.axis('equal')
12     
13 def throw_darts(r=1):
14     x = np.random.uniform(-r, r)
15     y = np.random.uniform(-r, r)
16     return (x, y)
17     
18 def plot_darts(x, y, color='b'):
19     plt.figure(1)
20     plt.plot(x, y, '.', mfc=color, mec=color)
21     plt.draw()
22     plt.pause(0.01)
23 
24 def plot_est(nums, ests):
25     plt.figure(2)
26     plt.plot(nums, ests, 'ko')
27     plt.plot(nums, m.pi*np.ones(len(nums)), 'g.')
28     plt.draw()
29     plt.pause(0.01)
30     
31 if __name__ == '__main__':
32     darts_to_throw = 10000
33     darts_thrown = 0
34     board_rad = 1
35     do_plot = 0
36     inside = 0
37     pi_est = []
38     x_i_vec = []
39     y_i_vec = []
40     x_o_vec = []
41     y_o_vec = []
42 
43     plt.figure(1)
44     plt.clf()
45     draw_board(r=board_rad)
46 
47     plt.figure(2)
48     plt.clf()
49 
50     while darts_thrown < darts_to_throw:
51         darts_thrown += 1
52         x_d, y_d = throw_darts(board_rad)
53         r_d = m.sqrt(x_d**2 + y_d**2)
54         color='m'
55         if r_d < board_rad:
56             inside += 1
57             x_i_vec += [x_d]
58             y_i_vec += [y_d]
59             color='g'
60         else:
61             x_o_vec += [x_d]
62             y_o_vec += [y_d]
63 
64         pi_est += [4*inside/darts_thrown]
65         print('{:5d}'.format(darts_thrown) + 
66                ' {:5d}'.format(inside) +
67                ' {:0.3e}'.format(pi_est[-1]))
68         
69         if do_plot:
70             plot_darts(x_d, y_d, color)
71             plot_est(darts_thrown, pi_est[-1])
72 
73     if not do_plot:
74         plot_darts(x_i_vec, y_i_vec, 'g')
75         plot_darts(x_o_vec, y_o_vec, 'm')
76         plot_est(np.linspace(1, darts_thrown, darts_thrown), pi_est)
77 
78     plt.figure(1)
79     draw_board(board_rad)

• Finally, you may only want to print the table every once in a while; you could set a variable called print_every equal to some value and then change the table printing code to something like:

        if (darts_thrown in [1, darts_to_throw] or 
            darts_thrown % print_every == 0):
            print('{:5d}'.format(darts_thrown) + 
                   ' {:5d}'.format(inside) +
                   ' {:0.3e}'.format(pi_est[-1]))

Note that print statement also needed to be indented to be under new if statement