In [1]:
Ntotal = 100
print(Ntotal)
In [51]:
Nin = 55.0
print(Ntotal/Nin)
In [52]:
import random
random.uniform(0.0,1.0)
Out[52]:
In [53]:
random.uniform(0.0,1.0)
Out[53]:
In [54]:
random.uniform(0.0,1.0)
Out[54]:
In [55]:
random.uniform(0.0,1.0)
Out[55]:
In [56]:
import math
import random
Ntotal = 100
Nin = 0
x = random.uniform(0.0,1.0)
y = random.uniform(0.0,1.0)
r = math.sqrt(x**2 + y**2)
print(x)
print(y)
print(r)
If you want to generate a series of numbers in a sequence, e.g. 1,2,3,4..., then you want the range() immutable sequence type (or its equivalent in some other Python library). See the example below. range() doesn't directly create a list of numbers; it is its own type in Python. To get a list from it, see below.
In [57]:
list_of_numbers=list(range(1,5))
print(list_of_numbers)
# Note that the list EXCLUDES the endpoint. If you want to get a list of 5 numbers, from 1-5,
# then you need to extend the endpoint by 1:
list_of_numbers=list(range(1,6))
print(list_of_numbers)
# or this
list_of_numbers=list(range(0,5))
print(list_of_numbers)
Let's put things together now. We can create a variable that stores the number of iterations ("loops") of the calculation we want to do. Let's call that Ntotal. Let's then loop 100 times and each time make a random x, random y, and from that compute $r=\sqrt{x^2+y^2}$.
Note that Python uses indentation to indicate a related block of code acting as a "subroutine" - a program within the program.
In [58]:
import math
Ntotal = 100
Nin = 0
for i in range(0,Ntotal):
x = random.uniform(0.0,1.0)
y = random.uniform(0.0,1.0)
r = math.sqrt(x**2 + y**2)
print("x=%f, y=%f, r=%f" % (x,y,r))
Now that we can make random "dots" in x,y and compute the radius (relative to 0,0) of each dot, let's employ "Accept/Reject" to see if something is in/on the circle or outside the circle. All we have to do is, for each $r$, test whether $r \le R$ or $r > R$. If the former, we have a dot in the circle - a hit! If the latter, then we have a dot out of the circle - a miss! We accept the hits and reject the misses. The total number of points define all the moves in the game, and the ratio of hits to the total will tell us about the area of the circle, and thus get us closer to $\pi$.
In [59]:
Ntotal = 100
Nin = 0
R = 1.0
for i in range(0,Ntotal):
x = random.uniform(0.0,1.0)
y = random.uniform(0.0,1.0)
r = math.sqrt(x**2 + y**2)
if r <= R:
Nin = Nin + 1
# alternatively, Nin += 1 (auto-increment by 1)
print("Number of dots: %d" % Ntotal)
print("Number of hits: %d" % Nin)
print("Number of misses: %d" % (Ntotal-Nin))
my_pi = 4.0*float(Nin)/float(Ntotal)
print("pi = %f" % (my_pi))
The number above is probably close to what you know as $\pi$, but likely not very precise. After all, we only threw 100 dots. We can increase precision by increasing the number of dots. In the code block below, feel free to play with Ntotal, trying different values. Observe how the computed value of $\pi$ changes. Would you say it's "closing in" on the value you know, or not? Try Ntotal at 1000, 5000, 10000, 50000, and 100000.
In the code block below, I have also added a computation of statistical error. The error should be a binomial error. Binomial errors occur when you have a bunch of things and you can classify them in two ways: as $A$ and $\bar{A}$ ("A" and "Not A"). In our case, they are hits or "not hits" (misses). So binomial errors apply. A binomial error is computed below.
In [60]:
Ntotal = 100
Nin = 0
R = 1.0
for i in range(0,Ntotal):
x = random.uniform(0.0,1.0)
y = random.uniform(0.0,1.0)
r = math.sqrt(x**2 + y**2)
if r <= R:
Nin = Nin + 1
# alternatively, Nin += 1 (auto-increment by 1)
my_pi = 4.0*float(Nin)/float(Ntotal)
my_pi_uncertainty = my_pi * math.sqrt(1.0/float(Nin) + 1.0/float(Ntotal))
print("pi = %.6f +/- %.6f (percent error= %.2f%%)" % (my_pi, my_pi_uncertainty, 100.0*my_pi_uncertainty/my_pi))