In [1]:
#%autosave 0
from IPython.core.display import HTML, display
display(HTML("<style>.container { width:100% !important; }</style>"))


Computing $\pi$ with the Monte-Carlo-Method


In [2]:
import random as rnd
import math

The unit circle $U$ is defined as the set $$U := \bigl\{ (x,y) \in \mathbb{R}^2 \bigm| x^2 + y^2 \leq 1 \bigr\}.$$ The set $U$ contains those points $(x,y)$ that have distance of $1$ or less from the origin $(0,0)$. The unit circle is a subset of the square $Q$ that is defined as $$Q := \bigl\{ (x,y) \in \mathbb{R}^2 \bigm| -1 \leq x \leq +1 \wedge -1 \leq y \leq +1 \bigr\}.$$ A simple algorithm to compute $\pi$ works as follows: We randomly create $n$ points $(x,y) \in Q$. Then we count the number of points that end up in the unit circle $U$. Call this number $k$. It is reasonable to assume that approximately $k$ is to $n$ as the area of $U$ is to the area of $Q$. As the area of $Q$ is $2 \cdot 2$ and the area of $U$ equals $\pi \cdot 1^2$, we should have $$\frac{k}{n} \approx \frac{\pi}{4}.$$ Multiplying by $4$ we get $$\pi \approx 4 \cdot \frac{k}{n}.$$ The function $\texttt{approximate_pi}(n)$ creates $n$ random points in $Q$ and approximates $\pi$ as $4 \cdot \frac{k}{n}$.


In [3]:
def approximate_pi(n):
    k = 0
    for _ in range(n):
        x = 2 * rnd.random() - 1
        y = 2 * rnd.random() - 1
        r = x * x + y * y
        if r <= 1:
            k += 1
    return 4 * k / n

Given a list $L = [x_1, \cdots, x_n]$, the function $\texttt{std_and_mean}(L)$ computes the pair $(\mu, \sigma)$, where $\sigma$ is the sample standard deviation of $L$, while $\mu$ is the mean of $L$. $\mu$ and $\sigma$ are defined as follows: $$ \mu = \sum\limits_{i=1}^n x_i$$ $$ \sigma = \sqrt{\sum\limits_{i=1}^n \frac{(x_i - \mu)^2}{N-1}} $$


In [4]:
def std_and_mean(L):
    N    = len(L)
    mean = sum(L) / N
    ss    = 0
    for x in L:
        ss += (x - mean) ** 2 
    ss /= (N - 1)
    std = math.sqrt(ss)
    return mean, std

The method $\texttt{confidence_interval}(k, n)$ runs $k$ approximations of $\pi$ using $n$ trials in each approximation run. It computes a $97.3\%$ confidence interval for the value of $\pi$.


In [5]:
def confidence_interval(k, n):
    L = []
    for _ in range(k):
        L.append(approximate_pi(n))
    𝜇, 𝜎 = std_and_mean(L)
    return 𝜇 - 3 * 𝜎, 𝜇, 𝜇 + 3 * 𝜎

In [6]:
%%time
n = 100
while n <= 10000000:
    lower, pi, upper = confidence_interval(100, n)
    print('%9d: %6f < 𝜋 < %6f, 𝜋 ≈ %6f, error: %6f' % (n, lower, upper, pi, abs(pi - math.pi)))
    n *= 10


      100: 2.630493 < 𝜋 < 3.645507, 𝜋 ≈ 3.138000, error: 0.003593
     1000: 2.987113 < 𝜋 < 3.297847, 𝜋 ≈ 3.142480, error: 0.000887
    10000: 3.092212 < 𝜋 < 3.190484, 𝜋 ≈ 3.141348, error: 0.000245
   100000: 3.124077 < 𝜋 < 3.156233, 𝜋 ≈ 3.140155, error: 0.001438
  1000000: 3.137091 < 𝜋 < 3.145157, 𝜋 ≈ 3.141124, error: 0.000469
 10000000: 3.140086 < 𝜋 < 3.143055, 𝜋 ≈ 3.141571, error: 0.000022
CPU times: user 6min 37s, sys: 469 ms, total: 6min 37s
Wall time: 6min 38s

In [ ]: