Probability concepts using Python

Dr. Tirthajyoti Sarkar, Fremont, CA 94536


This notebook illustrates the concept of probability (frequntist definition) using simple scripts and functions.

Set theory basics

Set theory is a branch of mathematical logic that studies sets, which informally are collections of objects. Although any type of object can be collected into a set, set theory is applied most often to objects that are relevant to mathematics. The language of set theory can be used in the definitions of nearly all mathematical objects.

Set theory is commonly employed as a foundational system for modern mathematics.

Python offers a native data structure called set, which can be used as a proxy for a mathematical set for almost all purposes.

In [1]:
# Directly with curly braces
Set1 = {1,2}
print (Set1)
print(type(Set1))
{1, 2}
<class 'set'>
In [2]:
my_list=[1,2,3,4]
my_set_from_list = set(my_list)
print(my_set_from_list)
{1, 2, 3, 4}

Membership testing with in and not in

In [3]:
my_set = set([1,3,5])
print("Here is my set:",my_set)
print("1 is in the set:",1 in my_set)
print("2 is in the set:",2 in my_set)
print("4 is NOT in the set:",4 not in my_set)
Here is my set: {1, 3, 5}
1 is in the set: True
2 is in the set: False
4 is NOT in the set: True

Set relations

  • Subset
  • Superset
  • Disjoint
  • Universal set
  • Null set
In [4]:
Univ = set([x for x in range(11)])
Super = set([x for x in range(11) if x%2==0])
disj = set([x for x in range(11) if x%2==1])
Sub = set([4,6])
Null = set([x for x in range(11) if x>10])
In [5]:
print("Universal set (all the positive integers up to 10):",Univ)
print("All the even positive integers up to 10:",Super)
print("All the odd positive integers up to 10:",disj)
print("Set of 2 elements, 4 and 6:",Sub)
print("A null set:", Null)
Universal set (all the positive integers up to 10): {0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10}
All the even positive integers up to 10: {0, 2, 4, 6, 8, 10}
All the odd positive integers up to 10: {1, 3, 5, 7, 9}
Set of 2 elements, 4 and 6: {4, 6}
A null set: set()
In [6]:
print('Is "Super" a superset of "Sub"?',Super.issuperset(Sub))
print('Is "Super" a subset of "Univ"?',Super.issubset(Univ))
print('Is "Sub" a superset of "Super"?',Sub.issuperset(Super))
print('Is "Super" disjoint with "disj"?',Sub.isdisjoint(disj))
Is "Super" a superset of "Sub"? True
Is "Super" a subset of "Univ"? True
Is "Sub" a superset of "Super"? False
Is "Super" disjoint with "disj"? True

Set algebra/Operations

  • Equality
  • Intersection
  • Union
  • Complement
  • Difference
  • Cartesian product
In [7]:
S1 = {1,2}
S2 = {2,2,1,1,2}
print ("S1 and S2 are equal because order or repetition of elements do not matter for sets\nS1==S2:", S1==S2)
S1 and S2 are equal because order or repetition of elements do not matter for sets
S1==S2: True
In [8]:
S1 = {1,2,3,4,5,6}
S2 = {1,2,3,4,0,6}
print ("S1 and S2 are NOT equal because at least one element is different\nS1==S2:", S1==S2)
S1 and S2 are NOT equal because at least one element is different
S1==S2: False

In mathematics, the intersection A ∩ B of two sets A and B is the set that contains all elements of A that also belong to B (or equivalently, all elements of B that also belong to A), but no other elements. Formally,

$$ {\displaystyle A\cap B=\{x:x\in A{\text{ and }}x\in B\}.} $$

3 sets intersection

In [9]:
# Define a set using list comprehension
S1 = set([x for x in range(1,11) if x%3==0])
print("S1:", S1)
S1: {9, 3, 6}
In [10]:
S2 = set([x for x in range(1,7)])
print("S2:", S2)
S2: {1, 2, 3, 4, 5, 6}
In [11]:
# Both intersection method or & can be used
S_intersection = S1.intersection(S2)
print("Intersection of S1 and S2:", S_intersection)

S_intersection = S1 & S2
print("Intersection of S1 and S2:", S_intersection)
Intersection of S1 and S2: {3, 6}
Intersection of S1 and S2: {3, 6}
In [12]:
S3 = set([x for x in range(6,10)])
print("S3:", S3)
S1_S2_S3 = S1.intersection(S2).intersection(S3)
print("Intersection of S1, S2, and S3:", S1_S2_S3)
S3: {8, 9, 6, 7}
Intersection of S1, S2, and S3: {6}

In set theory, the union (denoted by ∪) of a collection of sets is the set of all elements in the collection. It is one of the fundamental operations through which sets can be combined and related to each other. Formally,

$$ {A\cup B=\{x:x\in A{\text{ or }}x\in B\}} $$

Union of sets

In [13]:
# Both union method or | can be used
S1 = set([x for x in range(1,11) if x%3==0])
print("S1:", S1)
S2 = set([x for x in range(1,5)])
print("S2:", S2)

S_union = S1.union(S2)
print("Union of S1 and S2:", S_union)
S_union = S1 | S2
print("Union of S1 and S2:", S_union)
S1: {9, 3, 6}
S2: {1, 2, 3, 4}
Union of S1 and S2: {1, 2, 3, 4, 6, 9}
Union of S1 and S2: {1, 2, 3, 4, 6, 9}

Set algebra laws

Commutative law: $$ {\displaystyle A\cap B=B\cap A} $$ $$ {\displaystyle A\cup (B\cup C)=(A\cup B)\cup C} $$

Associative law: $$ {\displaystyle (A\cap B)\cap C=A\cap (B\cap C)} $$ $$ {\displaystyle A\cap (B\cup C)=(A\cap B)\cup (A\cap C)} $$

Distributive law: $$ {\displaystyle A\cap (B\cup C)=(A\cap B)\cup (A\cap C)} $$ $$ {\displaystyle A\cup (B\cap C)=(A\cup B)\cap (A\cup C)} $$

Complement

If A is a set, then the absolute complement of A (or simply the complement of A) is the set of elements not in A. In other words, if U is the universe that contains all the elements under study, and there is no need to mention it because it is obvious and unique, then the absolute complement of A is the relative complement of A in U. Formally,

$$ {\displaystyle A^{\complement }=\{x\in U\mid x\notin A\}.} $$

You can take the union of two sets and if that is equal to the universal set (in the context of your problem), then you have found the right complement

In [14]:
S=set([x for x in range (21) if x%2==0])
print ("S is the set of even numbers between 0 and 20:", S)
S is the set of even numbers between 0 and 20: {0, 2, 4, 6, 8, 10, 12, 14, 16, 18, 20}
In [15]:
S_complement = set([x for x in range (21) if x%2!=0])
print ("S_complement is the set of odd numbers between 0 and 20:", S_complement)
S_complement is the set of odd numbers between 0 and 20: {1, 3, 5, 7, 9, 11, 13, 15, 17, 19}
In [16]:
print ("Is the union of S and S_complement equal to all numbers between 0 and 20?", 
       S.union(S_complement)==set([x for x in range (21)]))
Is the union of S and S_complement equal to all numbers between 0 and 20? True

De Morgan's laws

$$ {\displaystyle \left(A\cup B\right)^{\complement }=A^{\complement }\cap B^{\complement }.} $$$$ {\displaystyle \left(A\cap B\right)^{\complement }=A^{\complement }\cup B^{\complement }.} $$

Complement laws

$$ {\displaystyle A\cup A^{\complement }=U.} $$$$ {\displaystyle A\cap A^{\complement }=\varnothing .} $$$$ {\displaystyle \varnothing ^{\complement }=U.} $$$$ {\displaystyle U^{\complement }=\varnothing .} $$$$ {\displaystyle {\text{If }}A\subset B{\text{, then }}B^{\complement }\subset A^{\complement }.} $$

Difference between sets

If A and B are sets, then the relative complement of A in B, also termed the set-theoretic difference of B and A, is the set of elements in B but not in A.

$$ {\displaystyle B\setminus A=\{x\in B\mid x\notin A\}.} $$

Set difference

In [17]:
S1 = set([x for x in range(31) if x%3==0])
print ("Set S1:", S1)
Set S1: {0, 3, 6, 9, 12, 15, 18, 21, 24, 27, 30}
In [18]:
S2 = set([x for x in range(31) if x%5==0])
print ("Set S2:", S2)
Set S2: {0, 5, 10, 15, 20, 25, 30}
In [19]:
S_difference = S2-S1
print("Difference of S1 and S2 i.e. S2\S1:", S_difference)

S_difference = S1.difference(S2)
print("Difference of S2 and S1 i.e. S1\S2:", S_difference)
Difference of S1 and S2 i.e. S2\S1: {25, 10, 20, 5}
Difference of S2 and S1 i.e. S1\S2: {3, 6, 9, 12, 18, 21, 24, 27}

Following identities can be obtained with algebraic manipulation:

$$ {\displaystyle C\setminus (A\cap B)=(C\setminus A)\cup (C\setminus B)} $$$$ {\displaystyle C\setminus (A\cup B)=(C\setminus A)\cap (C\setminus B)} $$$$ {\displaystyle C\setminus (B\setminus A)=(C\cap A)\cup (C\setminus B)} $$$$ {\displaystyle C\setminus (C\setminus A)=(C\cap A)} $$$$ {\displaystyle (B\setminus A)\cap C=(B\cap C)\setminus A=B\cap (C\setminus A)} $$$$ {\displaystyle (B\setminus A)\cup C=(B\cup C)\setminus (A\setminus C)} $$
  
$$ {\displaystyle A\setminus A=\emptyset} $$$$ {\displaystyle \emptyset \setminus A=\emptyset } $$$$ {\displaystyle A\setminus \emptyset =A} $$$$ {\displaystyle A\setminus U=\emptyset } $$

Symmetric difference

In set theory, the symmetric difference, also known as the disjunctive union, of two sets is the set of elements which are in either of the sets and not in their intersection. $$ {\displaystyle A\,\triangle \,B=\{x:(x\in A)\oplus (x\in B)\}}$$

$$ {\displaystyle A\,\triangle \,B=(A\smallsetminus B)\cup (B\smallsetminus A)} $$$${\displaystyle A\,\triangle \,B=(A\cup B)\smallsetminus (A\cap B)} $$

Some properties,

$$ {\displaystyle A\,\triangle \,B=B\,\triangle \,A,} $$$$ {\displaystyle (A\,\triangle \,B)\,\triangle \,C=A\,\triangle \,(B\,\triangle \,C).} $$

The empty set is neutral, and every set is its own inverse:

$$ {\displaystyle A\,\triangle \,\varnothing =A,} $$$$ {\displaystyle A\,\triangle \,A=\varnothing .} $$
In [20]:
print("S1",S1)
print("S2",S2)
print("Symmetric difference", S1^S2)
print("Symmetric difference", S2.symmetric_difference(S1))
S1 {0, 3, 6, 9, 12, 15, 18, 21, 24, 27, 30}
S2 {0, 5, 10, 15, 20, 25, 30}
Symmetric difference {3, 5, 6, 9, 10, 12, 18, 20, 21, 24, 25, 27}
Symmetric difference {3, 5, 6, 9, 10, 12, 18, 20, 21, 24, 25, 27}

Cartesian product

In set theory (and, usually, in other parts of mathematics), a Cartesian product is a mathematical operation that returns a set (or product set or simply product) from multiple sets. That is, for sets A and B, the Cartesian product A × B is the set of all ordered pairs (a, b) where a ∈ A and b ∈ B.

$$ {\displaystyle A\times B=\{\,(a,b)\mid a\in A\ {\mbox{ and }}\ b\in B\,\}.} $$

More generally, a Cartesian product of n sets, also known as an n-fold Cartesian product, can be represented by an array of n dimensions, where each element is an n-tuple. An ordered pair is a 2-tuple or couple. The Cartesian product is named after René Descartes whose formulation of analytic geometry gave rise to the concept.

In [21]:
A = set(['a','b','c'])
S = {1,2,3}
In [22]:
def cartesian_product(S1,S2):
    result = set()
    for i in S1:
        for j in S2:
            result.add(tuple([i,j]))
    return (result)
In [23]:
C = cartesian_product(A,S)
print("Cartesian product of A and S\n{} X {}:{}".format(A,S,C))
Cartesian product of A and S
{'b', 'c', 'a'} X {1, 2, 3}:{('a', 1), ('a', 2), ('b', 1), ('a', 3), ('b', 2), ('b', 3), ('c', 1), ('c', 2), ('c', 3)}
In [24]:
print("Length of the Cartesian product set:",len(C))
Length of the Cartesian product set: 9

Note that because these are ordered pairs, same element can be repeated inside the pair i.e. even if two sets contain some identical elements, they can be paired up in the Cartesian product.

Instead of writing functions ourselves, we could use the itertools library of Python. Remember to turn the resulting product object into a list for viewing and subsequent processing.

In [25]:
from itertools import product as prod
In [26]:
A = set([x for x in range(1,7)])
B = set([x for x in range(1,7)])
p=list(prod(A,B))

print("A is set of all possible throws of a dice:",A)
print("B is set of all possible throws of a dice:",B)
print ("\nProduct of A and B is the all possible combinations of A and B thrown together:\n",p)
A is set of all possible throws of a dice: {1, 2, 3, 4, 5, 6}
B is set of all possible throws of a dice: {1, 2, 3, 4, 5, 6}

Product of A and B is the all possible combinations of A and B thrown together:
 [(1, 1), (1, 2), (1, 3), (1, 4), (1, 5), (1, 6), (2, 1), (2, 2), (2, 3), (2, 4), (2, 5), (2, 6), (3, 1), (3, 2), (3, 3), (3, 4), (3, 5), (3, 6), (4, 1), (4, 2), (4, 3), (4, 4), (4, 5), (4, 6), (5, 1), (5, 2), (5, 3), (5, 4), (5, 5), (5, 6), (6, 1), (6, 2), (6, 3), (6, 4), (6, 5), (6, 6)]

Cartesian Power

The Cartesian square (or binary Cartesian product) of a set X is the Cartesian product $X^2 = X × X$. An example is the 2-dimensional plane $R^2 = R × R$ where R is the set of real numbers: $R^2$ is the set of all points (x,y) where x and y are real numbers (see the Cartesian coordinate system).

The cartesian power of a set X can be defined as:

${\displaystyle X^{n}=\underbrace {X\times X\times \cdots \times X} _{n}=\{(x_{1},\ldots ,x_{n})\ |\ x_{i}\in X{\text{ for all }}i=1,\ldots ,n\}.} $

The cardinality of a set is the number of elements of the set. Cardinality of a Cartesian power set is $|S|^{n}$ where |S| is the cardinality of the set S and n is the power.

We can easily use itertools again for calculating Cartesian power. The repeat parameter is used as power.

In [27]:
A = {'Head','Tail'} # 2 element set
p2=list(prod(A,repeat=2)) # Power set of power 2
print("Cartesian power 2 with length {}: {}".format(len(p2),p2))
print()
p3=list(prod(A,repeat=3)) # Power set of power 3
print("Cartesian power 3 with length {}: {}".format(len(p3),p3))
Cartesian power 2 with length 4: [('Tail', 'Tail'), ('Tail', 'Head'), ('Head', 'Tail'), ('Head', 'Head')]

Cartesian power 3 with length 8: [('Tail', 'Tail', 'Tail'), ('Tail', 'Tail', 'Head'), ('Tail', 'Head', 'Tail'), ('Tail', 'Head', 'Head'), ('Head', 'Tail', 'Tail'), ('Head', 'Tail', 'Head'), ('Head', 'Head', 'Tail'), ('Head', 'Head', 'Head')]

Permutations

In mathematics, the notion of permutation relates to the act of arranging all the members of a set into some sequence or order, or if the set is already ordered, rearranging (reordering) its elements, a process called permuting. The study of permutations of finite sets is a topic in the field of combinatorics.

We find the number of $k$-permutations of $A$, first by determining the set of permutations and then by calculating $\frac{|A|!}{(|A|-k)!}$. We first consider the special case of $k=|A|$, which is equivalent to finding the number of ways of ordering the elements of $A$.

In [28]:
import itertools
In [29]:
A = {'Red','Green','Blue'}
In [30]:
# Find all permutations of A
permute_all = set(itertools.permutations(A))
print("Permutations of {}".format(A))
print("-"*50)
for i in permute_all:
    print(i)
print("-"*50)
print;print ("Number of permutations: ", len(permute_all))
Permutations of {'Blue', 'Red', 'Green'}
--------------------------------------------------
('Red', 'Green', 'Blue')
('Green', 'Red', 'Blue')
('Blue', 'Red', 'Green')
('Green', 'Blue', 'Red')
('Red', 'Blue', 'Green')
('Blue', 'Green', 'Red')
--------------------------------------------------
Number of permutations:  6
In [31]:
import matplotlib.patches as mpatches
import matplotlib.pyplot as plt
In [32]:
from math import factorial
In [33]:
print("Factorial of 3:", factorial(3))
Factorial of 3: 6

Selecting k items out of a set containing n items and permuting

In [34]:
A = {'Red','Green','Blue','Violet'}
k=2
In [35]:
n = len(A)
permute_k = list(itertools.permutations(A, k))
print("{}-permutations of {}:  ".format(k,A))
print("-"*50)
for i in permute_k:
    print(i)
print("-"*50)
print ("Size of the permutation set = {}!/({}-{})! = {}".format(n,n,k, len(permute_k)))
2-permutations of {'Blue', 'Red', 'Green', 'Violet'}:  
--------------------------------------------------
('Blue', 'Red')
('Blue', 'Green')
('Blue', 'Violet')
('Red', 'Blue')
('Red', 'Green')
('Red', 'Violet')
('Green', 'Blue')
('Green', 'Red')
('Green', 'Violet')
('Violet', 'Blue')
('Violet', 'Red')
('Violet', 'Green')
--------------------------------------------------
Size of the permutation set = 4!/(4-2)! = 12
In [36]:
factorial(4)/(factorial(4-2))
Out[36]:
12.0

Combinations

Combinatorics is an area of mathematics primarily concerned with counting, both as a means and an end in obtaining results, and certain properties of finite structures. It is closely related to many other areas of mathematics and has many applications ranging from logic to statistical physics, from evolutionary biology to computer science, etc.

Combinatorics is well known for the breadth of the problems it tackles. Combinatorial problems arise in many areas of pure mathematics, notably in algebra, probability theory, topology, and geometry, as well as in its many application areas. Many combinatorial questions have historically been considered in isolation, giving an ad hoc solution to a problem arising in some mathematical context. In the later twentieth century, however, powerful and general theoretical methods were developed, making combinatorics into an independent branch of mathematics in its own right. One of the oldest and most accessible parts of combinatorics is graph theory, which by itself has numerous natural connections to other areas. Combinatorics is used frequently in computer science to obtain formulas and estimates in the analysis of algorithms.

We find the number of $k$-combinations of $A$, first by determining the set of combinations and then by simply calculating: $$\frac{|A|!}{k!\times(|A|-k)!}$$

In combination, order matters, unlike permutations.

In [37]:
# Print all the k-combinations of A
choose_k = list(itertools.combinations(A,k))
print("%i-combinations of %s:  " %(k,A))
for i in choose_k:
    print(i)
print;print("Number of combinations = %i!/(%i!(%i-%i)!) = %i" %(n,k,n,k,len(choose_k)  ))
2-combinations of {'Blue', 'Red', 'Green', 'Violet'}:  
('Blue', 'Red')
('Blue', 'Green')
('Blue', 'Violet')
('Red', 'Green')
('Red', 'Violet')
('Green', 'Violet')
Number of combinations = 4!/(2!(4-2)!) = 6

Putting it all together - some probability calculation examples

Problem 1: Two dice

Two dice are rolled together. What is the probability of getting a total which is a multiple of 3?

In [38]:
n_dice = 2
dice_faces = {1,2,3,4,5,6}

# Construct the event space i.e. set of ALL POSSIBLE events
event_space = set(prod(dice_faces,repeat=n_dice))
In [39]:
for outcome in event_space:
    print(outcome,end=', ')
(1, 3), (6, 6), (5, 6), (2, 1), (6, 2), (1, 6), (5, 1), (2, 5), (1, 2), (3, 3), (5, 5), (4, 4), (6, 3), (1, 5), (3, 6), (2, 2), (4, 1), (1, 1), (6, 4), (3, 2), (2, 6), (5, 4), (4, 5), (5, 2), (1, 4), (2, 3), (4, 2), (6, 5), (3, 5), (5, 3), (4, 6), (6, 1), (3, 1), (4, 3), (3, 4), (2, 4), 
In [40]:
# What is the set we are interested in?
favorable_outcome = []
for outcome in event_space:
    x,y = outcome
    if (x+y)%3==0:
        favorable_outcome.append(outcome)
favorable_outcome = set(favorable_outcome)
In [41]:
for f_outcome in favorable_outcome:
    print(f_outcome,end=', ')
(1, 2), (5, 4), (3, 3), (6, 6), (4, 5), (2, 1), (6, 3), (1, 5), (3, 6), (4, 2), (5, 1), (2, 4), 
In [42]:
prob = len(favorable_outcome)/len(event_space)
In [43]:
print("The probability of getting a sum which is a multiple of 3 is: ", prob)
The probability of getting a sum which is a multiple of 3 is:  0.3333333333333333

Problem 2: Five dice!

Five dice are rolled together. What is the probability of getting a total which is a multiple of 5 but not a multiple of 3?

In [44]:
n_dice = 5
dice_faces = {1,2,3,4,5,6}

# Construct the event space i.e. set of ALL POSSIBLE events
event_space = set(prod(dice_faces,repeat=n_dice))
In [45]:
6**5
Out[45]:
7776
In [46]:
# What is the set we are interested in?
favorable_outcome = []
for outcome in event_space:
    d1,d2,d3,d4,d5 = outcome
    if (d1+d2+d3+d4+d5)%5==0 and (d1+d2+d3+d4+d5)%3!=0 :
        favorable_outcome.append(outcome)
favorable_outcome = set(favorable_outcome)
In [47]:
prob = len(favorable_outcome)/len(event_space)
print("The probability of getting a sum, which is a multiple of 5 but not a multiple of 3, is: ", prob)
The probability of getting a sum, which is a multiple of 5 but not a multiple of 3, is:  0.11625514403292181

Problem 2 solved using set difference

In [48]:
multiple_of_5 = []
multiple_of_3 = []

for outcome in event_space:
    d1,d2,d3,d4,d5 = outcome
    if (d1+d2+d3+d4+d5)%5==0:
        multiple_of_5.append(outcome)
    if (d1+d2+d3+d4+d5)%3==0:
        multiple_of_3.append(outcome)

favorable_outcome = set(multiple_of_5).difference(set(multiple_of_3))
In [49]:
for i in list(favorable_outcome)[:5]:
    a1,a2,a3,a4,a5=i
    print("{}, SUM: {}".format(i,a1+a2+a3+a4+a5))    
(3, 1, 1, 1, 4), SUM: 10
(4, 1, 6, 3, 6), SUM: 20
(6, 2, 2, 4, 6), SUM: 20
(1, 4, 6, 6, 3), SUM: 20
(5, 6, 4, 5, 5), SUM: 25
In [50]:
prob = len(favorable_outcome)/len(event_space)
print("The probability of getting a sum, which is a multiple of 5 but not a multiple of 3, is: ", prob)
The probability of getting a sum, which is a multiple of 5 but not a multiple of 3, is:  0.11625514403292181

Computing pi ($\pi$) with a random dart throwing game and using probability concept

The number $\pi$ is a mathematical constant. Originally defined as the ratio of a circle's circumference to its diameter, it now has various equivalent definitions and appears in many formulas in all areas of mathematics and physics. It is approximately equal to 3.14159. It has been represented by the Greek letter "$\pi$" since the mid-18th century, though it is also sometimes spelled out as "pi". It is also called Archimedes' constant.

Being an irrational number, $\pi$ cannot be expressed as a common fraction (equivalently, its decimal representation never ends and never settles into a permanently repeating pattern).

What is the logic behind computing $\pi$ by throwing dart randomly?

Imagine a square dartboard.

Then, the dartboard with a circle drawn inside it touching all its sides.

And then, you throw darts at it. Randomly. That means some fall inside the circle, some outside. But assume that no dart falls outside the board.

boards

At the end of your dart throwing session,

  • You count the fraction of darts that fell inside the circle of the total number of darts thrown.
  • Multiply that number by 4.
  • The resulting number should be pi. Or, a close approximation if you had thrown a lot of darts.

The idea is extremely simple. If you throw a large number of darts, then the probability of a dart falling inside the circle is just the ratio of the area of the circle to that of the area of the square board. With the help of basic mathematics, you can show that this ratio turns out to be $\frac{\pi}{4}$. So, to get $\pi$, you just multiply that number by 4.

The key here is to simulate the throwing of a lot of darts so as to make the fraction equal to the probability, an assertion valid only in the limit of a large number of trials of this random event. This comes from the law of large number or the frequentist definition of probability.

See also the concept of Buffon's Needle

In [51]:
from math import pi,sqrt
import random
import matplotlib.pyplot as plt
import numpy as np

Center point and the side of the square

In [52]:
# Center point
x,y = 0,0
# Side of the square
a = 2

Function to simulate a random throw of a dart aiming at the square

In [53]:
def throw_dart():
    """
    Simulates the randon throw of a dirt. It can land anywhere in the square (uniformly randomly)
    """
    # Random final landing position of the dirt between -a/2 and +a/2 around the center point
    position_x = x+a/2*(-1+2*random.random())
    position_y = y+a/2*(-1+2*random.random())
    
    return (position_x,position_y)
In [54]:
throw_dart()
Out[54]:
(0.22132208501841943, 0.7844172623010937)

Function to determine if the dart landed inside the circle

In [55]:
def is_within_circle(x,y):
    """
    Given the landing coordinate of a dirt, determines if it fell inside the circle
    """
    # Side of the square
    a = 2
    
    distance_from_center = sqrt(x**2+y**2)
    
    if distance_from_center < a/2:
        return True
    else:
        return False
In [56]:
is_within_circle(1.9,1.9)
Out[56]:
False
In [57]:
is_within_circle(0.2,-0.6)
Out[57]:
True

Now, throw a few darts

In [58]:
n_throws = 10
count_inside_circle=0
In [59]:
for i in range(n_throws):
    r1,r2=throw_dart()
    if is_within_circle(r1,r2):
        count_inside_circle+=1

Compute the ratio of count_inside_circle and n_throws

In [60]:
ratio = count_inside_circle/n_throws

Is it approximately equal to $\frac{\pi}{4}$?

In [61]:
print(4*ratio)
2.8

Not exactly. Let's try with a lot more darts!

In [62]:
n_throws = 10000
count_inside_circle=0

for i in range(n_throws):
    r1,r2=throw_dart()
    if is_within_circle(r1,r2):
        count_inside_circle+=1
In [63]:
ratio = count_inside_circle/n_throws
print(4*ratio)
3.136

Let's functionalize this process and run a number of times

In [64]:
def compute_pi_throwing_dart(n_throws):
    """
    Computes pi by throwing a bunch of darts at the square
    """
    n_throws = n_throws
    count_inside_circle=0
    for i in range(n_throws):
        r1,r2=throw_dart()
        if is_within_circle(r1,r2):
            count_inside_circle+=1
            
    result = 4*(count_inside_circle/n_throws)
    
    return result

Now let us run this experiment a few times and see what happens.

In [65]:
n_exp=[]
pi_exp=[]
n = [int(10**(0.5*i)) for i in range(1,15)]
for i in n:
    p = compute_pi_throwing_dart(i)
    pi_exp.append(p)
    n_exp.append(i)
    print("Computed value of pi by throwing {} darts is: {}".format(i,p))
Computed value of pi by throwing 3 darts is: 4.0
Computed value of pi by throwing 10 darts is: 2.4
Computed value of pi by throwing 31 darts is: 3.225806451612903
Computed value of pi by throwing 100 darts is: 3.44
Computed value of pi by throwing 316 darts is: 3.367088607594937
Computed value of pi by throwing 1000 darts is: 3.152
Computed value of pi by throwing 3162 darts is: 3.135989879822897
Computed value of pi by throwing 10000 darts is: 3.1004
Computed value of pi by throwing 31622 darts is: 3.141483777117197
Computed value of pi by throwing 100000 darts is: 3.13996
Computed value of pi by throwing 316227 darts is: 3.145550506439994
Computed value of pi by throwing 1000000 darts is: 3.142276
Computed value of pi by throwing 3162277 darts is: 3.140461129749228
Computed value of pi by throwing 10000000 darts is: 3.1421584
In [66]:
plt.figure(figsize=(8,5))
plt.title("Computing pi with \nincreasing number of random throws",fontsize=20)
plt.semilogx(n_exp, pi_exp,c='k',marker='o',lw=3)
plt.xticks(fontsize=14)
plt.yticks(fontsize=14)
plt.xlabel("Number of random throws",fontsize=15)
plt.ylabel("Computed value of pi",fontsize=15)
plt.hlines(y=3.14159,xmin=1,xmax=1e7,linestyle='--')
plt.text(x=10,y=3.05,s="Value of pi",fontsize=17)
plt.grid(True)
plt.show()