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