In [1]:

```
%load_ext watermark
```

In [2]:

```
watermark -d -v -a "Sebastian Raschka"
```

While I was listening to the recent episode of the Programming Throwdown podcast on my way home, I heard about this interesting concept of rejection sampling. This is one of those simple ideas combined with the powerful principles of statistics that I find quite fascinating.

At its core, rejection sampling is similar to the popular Monte Carlo sampling with the difference of an additional bound. The goal of rejection sampling is to simplify the task of drawing random samples from a complex probability distribution by using a uniform distribution instead; random samples drawn from the uniform distribution that lie outside certain boundary criteria are rejected, and all samples within the boundary are accepted, respectively.

A mathematical proof can be found here.

In [3]:

```
%matplotlib inline
import matplotlib.pyplot as plt
```

In [4]:

```
center = (4,4)
radius = 2
def plot_circle(center, radius):
""" Function to plot a circle. """
fig = plt.figure(figsize=(6,6))
circle = plt.Circle(center, radius, fill=False, color='b')
plt.ylim([0,8])
plt.xlim([0,8])
fgca = fig.gca()
fgca.add_artist(circle)
return fgca
plot_circle(center, radius)
plt.show()
```

Now, we can draw a simple square around the circle that will represent our uniform distribution.

In [5]:

```
from matplotlib.patches import Rectangle
def plot_square(center, radius):
""" Function to plot a square. """
fgca = plot_circle(center, radius)
fgca.add_patch(Rectangle((center[0] - radius, center[1] - radius),
2*radius, 2*radius, alpha=0.1))
return fgca
plot_square(center, radius)
plt.show()
```

In [6]:

```
import random
random.seed(567)
```

In [7]:

```
def gen_points(n, center, radius):
"""
Function that generates
n x,y coordinates in a square of length radius*2.
"""
x_coords = []
y_coords = []
for i in range(n):
x_coords.append(random.random()*center[0]+radius)
y_coords.append(random.random()*center[1]+radius)
return x_coords, y_coords
```

In [8]:

```
x, y = gen_points(1, center, radius)
fgca = plot_square(center, radius)
fgca.plot(x, y, linestyle="", marker="x", color="red")
plt.show()
```

Let us generate 1000 random points and check if our function works correctly.

In [9]:

```
x, y = gen_points(1000, center, radius)
fgca = plot_square(center, radius)
fgca.plot(x, y, linestyle="", marker="x", color="red")
plt.show()
```

The plot above looks fine. In the last step, we only need to reject those points lie outside the circle, which is pretty straight forward using a Euclidean distance measure, that is

\begin{equation} \sqrt{(x-y)^2} = |x-y|. \end{equation}In [10]:

```
def reject(radius, center, x_coords, y_coords):
""" Returns those coordinates that fall within the circle. """
x_clean = []
y_clean = []
for x,y in zip(x_coords, y_coords):
if ((x - center[0])**2 + (y-center[1])**2)**0.5 <= radius:
x_clean.append(x)
y_clean.append(y)
return x_clean, y_clean
```

In [11]:

```
x_clean, y_clean = reject(radius, center, x, y)
fgca = plot_square(center, radius)
fgca.plot(x_clean, y_clean, linestyle="", marker="x", color="red")
plt.show()
```

\begin{equation} A_{\text{rectangle}} = (2 \times R)^2 \end{equation}

\begin{equation} \hat{A}_{\text{est_circle}} = A_{\text{rectangle}} \times \frac{\text{# points inside circle}}{\text{# all points}} \end{equation}

In [12]:

```
def estimate_circle_area(n, center, radius):
""" Returns the estimated circle area via rejection sampling. """
rect_area = (2*radius)**2
x, y = gen_points(n, center, radius)
x_clean, y_clean = reject(radius, center, x, y)
est_circle_area = rect_area * (len(x_clean)/len(x))
return est_circle_area
```

In [13]:

```
print('Estimated circle area: %s' %estimate_circle_area(100000, center, radius))
```

Now, let's double-check how close we got using the more accurate equation:

\begin{equation} A_{\text{circle}} = \pi \times R^2 \end{equation}In [14]:

```
from math import pi
print('Circle area using pi: %s' %(pi*radius**2))
```

\begin{equation} \hat{\pi} = \frac{\hat{A}_{\text{est_circle}}}{R^2} \end{equation}

In [15]:

```
def approximate_pi(n, center, radius):
""" Returns an approximation of pi via rejection sampling. """
circ_area = estimate_circle_area(n, center, radius)
return circ_area/radius**2
for i in (10, 10**2, 10**4, 10**7):
pi_est = approximate_pi(i, center, radius)
print('Pi estimate: %s (n=%s)' %(pi_est, i))
```

Again, this estimate is pretty accurate up to the 4th digit after the decimal point.