Monte Carlo determination of $\pi$¶

Imagine that we have a circle of radius 1 ($r = 1$), centered at the point $(0, 0)$ in Cartesian coordinates.

We can further imagine that circle being enclosed by a square with the same center and with sides of length $2r$, in this case 2.

Given a square with sides of length $2r$, the square's area is

$$A_{square} = (2r)^2 = 4r^2$$

whereas the area of a circle with radius $r$ is $$A_{circle} = \pi r^2$$

Therefore the ratio of the area of the circle to that of the enclosing square is

$$\frac{A_{circle}}{A_{square}} = \frac{\pi r^2}{4r^2} = \frac{\pi}{4}$$

.

By rearranging this expression, we see we can define $\pi$ as

$$\pi = 4\frac{A_{circle}}{A_{square}}$$

And this suggests a way to calculate $\pi$!

If we're able to figure out the ratio of the areas of a circle and the smallest square that can hold that circle, we can figure out $\pi$. And fortunately, this ratio is equal to the probability that a randomly chosen point inside the square also falls inside the circle.

This means that we can calculate this ratio (and therefore $\pi$) using a Monte Carlo simulation by choosing a large number of random points inside a square and tracking how many fall inside the circle.

Pseudo-code¶

Given the above, our algorithm for determining $\pi$ looks like this:

For each of $N$ iterations,
1. Select a random point inside the square depicted above.
2. Determine if the point also falls inside the largest circle that can be enclosed within this square.
3. Keep track of whether or not this point fell inside the circle. At the end of $N$ iterations, you want to know $M$ -- the number of the $N$ random points that fell inside the circle!
Calculate $\pi$ as $4\frac{M}{N}$

Note: As you increase the value of $N$, you'll be able to estimate pi with greater accuracy!

Exercise¶

Write a function that calculates $\pi$ using Julia.

Feel free to get started immediately, or to make use of the practice exercises below which have you solve pieces of the problem one at a time.

In [35]:

function dist(x, y)
    return sqrt(x^2 + y^2)
end

Out[35]:

dist (generic function with 1 method)

In [36]:

function myrand()
   return rand([-1, 1]) * rand(), rand([-1, 1]) * rand() 
end

Out[36]:

myrand (generic function with 1 method)

In [37]:

function inside_circle(x, y)
    if dist(x, y) < 1.0
        return 1
    else
        return 0
    end
end

Out[37]:

inside_circle (generic function with 1 method)

In [38]:

"""
    calc_pi(N)
    
Calculate pi using a Monte Carlo simulation in which `N` sampling points are chosen within a square.
"""
function calc_pi(N)
    points_inside = 0
    for i in 1:N
        x, y = myrand()
        points_inside += inside_circle(x, y)
    end
    return 4 * (points_inside / N)
end

Out[38]:

calc_pi

Function practice¶

Exercise¶

Given variables x and y defining a point $(x, y)$, write a function dist(x, y) that returns the distance of $(x, y)$ from $(0, 0)$.

In [4]:

function dist(x, y)
    return sqrt(x^2 + y^2)
end

Out[4]:

dist (generic function with 1 method)

In [6]:

x, y = rand(), rand()

Out[6]:

(0.09177130543526113, 0.880104215610646)

In [7]:

dist(x, y)

Out[7]:

0.8848759251086687

Exercise¶

Write a function myrand that returns (x, y), where x and y are pseudo-random numbers uniformly distributed between -1 and +1.

Hint: The rand() function returns a random number between 0 and 1. Passing a collection as the first argument of rand will cause rand to return a random value from amongst the elements of the collection.

Ex. rand(animals) will return one of the strings contained in animals:

In [8]:

animals = ["cat", "dog", "bird"]
rand(animals)

Out[8]:

"cat"

In [9]:

function myrand()
   return rand([-1, 1]) * rand(), rand([-1, 1]) * rand() 
end

Out[9]:

myrand (generic function with 1 method)

In [13]:

myrand()

Out[13]:

(-0.9230968519773852, -0.38550538642272825)

Conditional practice¶

Exercise¶

Given variables x and y defining a point $(x, y)$, write a conditional statement that prints "Less than 1!" if the distance of $(x, y)$ from $(0, 0)$ is less than 1.

In [2]:

x, y = rand(2)

Out[2]:

2-element Array{Float64,1}:
 0.2753495990422792  
 0.011453969623016924

In [3]:

if sqrt(x^2 + y^2) < 1.0
    println("Less than 1!")
end

Less than 1!

In [5]:

if dist(x, y) < 1.0
    println("Less than 1!")
end

Less than 1!

Exercise¶

Write a function inside_circle(x, y) that takes coordinates/variables x and y and determines whether the point $(x, y)$ falls inside a circle of radius 1 that is centered at $(0, 0)$. Return 1 if the point falls within the circle and 0 if not.

Hint: To fall within the circle, the distance of the point to the circle's center should be less than the radius.

In [15]:

function inside_circle(x, y)
    if sqrt(x^2 + y^2) < 1.0
        return 1
    else
        return 0
    end
end

Out[15]:

inside_circle (generic function with 1 method)

In [14]:

inside_circle(x, y) = (sqrt(x^2 + y^2) < 1.0) ? 1 : 0

Out[14]:

inside_circle (generic function with 1 method)

Loop practice¶

Exercise¶

Use a loop to generate 100 random points $(x_i, y_i)$ where $-1.0 <= x <= 1.0$ and $-1.0 <= y <= 1.0$. Print the number of times the points fell within a circle of radius 1, centered at $(0, 0)$.

In [27]:

N = 0
for i in 1:100
    x, y = myrand()
    N += inside_circle(x, y)
end
println("$N points fell inside the circle!")

83 points fell inside the circle!

Without work that we've already done in previous exercises:

In [29]:

N = 0
for i in 1:100
    x = rand([-1, 1]) * rand()
    y = rand([-1, 1]) * rand()
    if sqrt(x^2 + y^2) < 1.0
        N += 1
    end
end
println("$N points fell inside the circle!")

78 points fell inside the circle!

Array practice¶

Exercise¶

Use an array comprehension to create an array of 100 points in 2D Cartesian space.

In [30]:

[myrand() for i in 1:100]

Out[30]:

100-element Array{Tuple{Float64,Float64},1}:
 (0.8235803117901142, -0.6736936601474774)  
 (0.3558640960890065, -0.7080181200055458)  
 (0.6161108313295864, 0.10128327704482909)  
 (0.06767258090743411, 0.48567762731410147) 
 (-0.01869616622547987, -0.4229478969051781)
 (0.8599279235169288, -0.20686834775240293) 
 (0.25769489684403424, -0.9587876421300074) 
 (-0.7108573976816956, -0.7344179861513818) 
 (0.7646953306782522, 0.33019531478167496)  
 (0.20513540415981146, -0.9106607051488109) 
 (-0.15357265925894148, -0.4578539987702732)
 (0.76617200964695, -0.7296079522154135)    
 (-0.5640083749041798, -0.8156093136312135) 
 ⋮                                          
 (0.4825477312267852, -0.6023767911672253)  
 (0.20483110741326693, 0.16210952105073861) 
 (-0.3645209544985353, -0.46411948118371793)
 (0.1762587416098942, -0.876202029020815)   
 (-0.4087320814383706, 0.3464881102351882)  
 (0.803515030181229, 0.227084405652221)     
 (-0.14559086789114306, -0.543519582178885) 
 (0.770377554719148, -0.5780172129728407)   
 (0.9773512106749624, 0.7160989375188618)   
 (0.5643735149375608, -0.9564993822734364)  
 (0.31253625667570106, 0.9430302711198242)  
 (-0.7897705751424646, -0.8060195845846594)

Exercise¶

Use an array comprehension to create an array of the distances of 100 points in 2D Cartesian space from a point at $(0, 0)$.

Hint: You may find "splat", i.e. ..., useful here. Adding ... after a collection passed as an input argument to a function "opens up" that collection, passing each of the elements of the collection as inputs to the function. For example, if

a = [1, 2, 3]
add_three_things(x, y, z) = x + y + z

then add_three_things(a) will give an error, but add_three_things(a...) will yield 6, just as if you had run add_three_things(a[1], a[2], a[3])

In [34]:

a = [1, 2, 3]
add_three_things(x, y, z) = x + y + z
add_three_things(a...)

Out[34]:

In [31]:

[dist(myrand()...) for i in 1:100]

Out[31]:

100-element Array{Float64,1}:
 0.9383872456685842 
 0.5839277168087241 
 0.7355502495040855 
 0.8441859785114507 
 0.6807536016582622 
 0.7632412130157962 
 0.7317234224262192 
 0.8824052199190965 
 0.7755182140496781 
 0.9744510181703258 
 0.25005959593047916
 1.0612231389451416 
 1.0606971366415012 
 ⋮                  
 0.9484156186801598 
 0.9270152199659667 
 0.6614786312111012 
 0.6478811880980543 
 0.23342156646680173
 1.1057076816888292 
 0.5958913950161429 
 1.2529131820213943 
 0.7256403533033896 
 0.910877586092377  
 0.779338060314894  
 0.16368015145789885

In [ ]: