These are some methods for uniformly random sampling in and on some common shapes including circles and spheres.

In [1]:
%matplotlib inline
import numpy as np
from numpy.random import uniform,standard_normal
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D


# Uniform sampling on a line¶

I will assume that we have access to a 1D uniform distribution sampling function and a 1D normal distribution sampling function. So, uniform sampling on a line becomes trivial by mapping the line to $[0,1]$. Then uniformly random draws from $[0,1]$ can be mapped back to the line.

# Uniform sampling in a Triangle¶

We plot a triangle to begin with and then generate $N=1000$ uniformly distributed random points inside it.

$P = (1 - \sqrt{r_1}) A + (\sqrt{r_1} (1 - r_2)) B + (r_2 \sqrt{r_1}) C$ where $r_1, r_2 \sim U[0, 1]$ gives us the required sampling.[src]

In [2]:
N=1000
verts=np.asarray([[0,0],[0.5,1],[1,0]])
fig,ax=plt.subplots(figsize=[6,6])

p1 = np.sqrt(uniform(0,1,N)) # taking sqrt here itself
p2 = uniform(0,1,N)

ts=np.empty([N,2])
ts[:,0] = (1-p1)*verts[0,0] + p1*(1-p2)*verts[1,0] + p2*p1*verts[2,0]
ts[:,1] = (1-p1)*verts[0,1] + p1*(1-p2)*verts[1,1] + p2*p1*verts[2,1]
_=ax.plot(ts[:,0],ts[:,1],'ro')


An alternative solution is to form a parallellogram using two sides and

# Uniform sampling in a Square¶

In [3]:
N=1000
verts=np.asarray([[0,0],[0,1],[1,1],[1,0]])
fig,ax=plt.subplots(figsize=[6,6])

u1 = uniform(0,1,N)
u2 = uniform(0,1,N)

_=ax.plot(u1,u2,'ro')


# Uniform sampling inside a circle¶

## Naive [WRONG] procedure¶

We want to generate $N$ random points inside the circle, such that they are distributed with an equal density. Using polar coordinates, $r$ and $\theta$, any point inside a circle or radius $r_0$ can be represented as

$$(x,y)=(r~cos(\theta),r~sin(\theta))$$ for $$r\in[0,r_0],~\theta\in[0,2\pi]$$

The naive approach is to randomly pick $N$ values each of $r\in[0,r_0]$ and a $\theta\in[0,2\pi]$ and calculate $x$ and $y$. As we shall see this does not give us a uniform distribution inside the circle. Lets say $N=1000$ and $r_0=2$.

In [4]:
r0 = 2.0
N  = 1000

r  = r0*uniform(0,1,N)
t  = uniform(0,2*np.pi,N)

x,y= r*np.cos(t),r*np.sin(t)

fig,ax=plt.subplots(figsize=(6,6))
ax.plot(x,y,'o')
ax.set_aspect('equal')


Clearly, we do not get a uniform distribution of points in the cirle this way. The points are clustered near the center.

## Best procedure¶

Take a step back and see why this went wrong. Lets say we already have a uniform sampling inside the circle, then on integrating over theta we would find the the density along $r$ is not constant (i.e not uniform). In fact, the density at any point would be proportional to the value of $r$ at that point. This means drawing $r$ from a uniform distribution was not right. We need to draw $r$ from a distribution that starts at 0 and increases linearly till $r=r_0$. Now that we know this, we can use the Inverse transform sampling procedure to draw from the correct distribution. For this to work, we need to integrate the $PDF$ to get the $CDF$ and then invert the $CDF$. Applying the inverse of the CDF on a uniformly sampled random variable will give us draws from the $PDF$ that we need. Given a continuous uniform variable $U$ in $[0,1]$ and an invertible cumulative distribution function $F_{X}$, the random variable $X=F_{X}^{-1}(U)$ has distribution $F_{X}$.[*].

In this particular case the $PDF=\frac{1}{r_0}\int{2r}dr$ and $CDF=r^2$ and so $CDF^{-1}(x)=x^{\frac{1}{2}}$. So, for the radius we draw uniform samples again but this time instead of using it directly, we take square-root. Alternative way to think about this is to imagine that we are taking uniform samples for $r^2$ from $[0,r_0^2]$.

In [5]:
r  = r0*np.sqrt(uniform(0,1,N))
t  = uniform(0,2*np.pi,N)
#################################
x,y= r*np.cos(t),r*np.sin(t)
fig,ax=plt.subplots(figsize=(6,6))
ax.plot(x,y,'o')
ax.set_aspect('equal')


This looks uniform now. This probably the best and most efficient way for sampling inside a circle. It needed $2N$ samples total.

## Alternative procedure¶

There is an alternative way to sample the "linear" probability distribution, without going the inverse transform sampling route. This needs two independent draws from a uniform distribution for the radius (i.e. three draws total). First, we add the two draws. The resulting points look like they have been drawn from a distribution that looks like a hat function with a peak at $r_0$ and domain from $0$ to $2r_0$. 'Folding back' the portion in $[r_0,2r_0]$ onto $[0,r_0]$, we get our linearly increasing distribution.

In [6]:
x1 = uniform(0.0,r0,N)
x2 = uniform(0.0,r0,N)
r  = x1+x2
r  = r0-np.abs(r-r0)# same as r=2*r0-r if r>r0 else r

########################################
x,y= r*np.cos(t),r*np.sin(t)
fig,ax=plt.subplots(figsize=(6,6))
ax.plot(x,y,'o')
ax.set_aspect('equal')


This looks alright as well. But, we needed to draw $3N$ samples.

## Rejection Sampling Procedure¶

This is a general alternative procedure to get uniform density in various circumstances. First, uniform samples in the square region $[-r_0,r_0]\times[-r_0,r_0]$ are drawn. Then, only the samples that fall inside the circle are kept. This is repeated until we have the necessary number of samples.

In [7]:
x=[]
y=[]
r0square=r0**2
ns=0
k=0
while(ns<N):
k+=1
xp = uniform(-r0,r0,1)
yp = uniform(-r0,r0,1)
if (xp**2+yp**2) < r0square:#acceptable
x.append(xp)
y.append(yp)
ns+=1

########################################
fig,ax=plt.subplots(figsize=(6,6))
ax.plot(x,y,'o')
ax.set_aspect('equal')
print('Total number of samples drawn:',2*k)

Total number of samples drawn: 2516


Unlike the previous methods, number of draws needed for this method is not known apriori. This is very simple to program easily extendable to 3D and different shapes like ellipses and ellipsoids for instance.

# Sampling on a sphere¶

Sampling 'on a circle' is trivial. We can draw from a uniform distribution in $[0,2\pi]$. This is the same as generating random directions or randomly oriented vectors in 2D world. The same in 3D is less clear. Generating a uniform sampling on the surface of a sphere is the same as generating uniformly randomly oriented vectors in 3D space. This is different from randomly rotating a 3d object. A vector does not change when rotated about itself, whereas a general 3d object would.

## Naive [WRONG] method¶

The surface of a sphere can be parametrized in $\theta$ and $\phi$. Lets see what we get when we naively sample these two from a uniform distribution. We'll take it to be a unit sphere.

In [8]:
N=500
phi   = np.random.uniform(0,2*np.pi,size=N)
theta = np.random.uniform(0,2*np.pi,size=N)

st=np.sin(theta)
ct=np.cos(theta)
sp=np.sin(phi)
cp=np.cos(phi)

v=np.c_[st*cp,st*sp,ct]

fig,ax=plt.subplots(figsize=(10,10))
ax = plt.axes(projection='3d')
ax.plot(v[:,0],v[:,1],v[:,2],'o')
ax.set_aspect('equal')

# Make spherical surface
u = np.linspace(0, 2 * np.pi, 100)
v = np.linspace(0, np.pi, 100)
x =  np.outer(np.cos(u), np.sin(v))
y =  np.outer(np.sin(u), np.sin(v))
z =  np.outer(np.ones(np.size(u)), np.cos(v))

# Plot the surface
collection=ax.plot_surface(x, y, z, color='r',alpha=0.5)


We see points clustering at the poles, similar to what was seen with the naive sampling of the circles. To fix this we can do inverse transform sampling of theta as shown below.

## Inverse transform Sampling¶

When a sphere is sliced along the latitudes(i.e. with parallel planes), the area of the slice does not depend on the latitude. This means that $z\in[-1,1]$ and $\phi\in[0,2\pi]$ can be picked from uniformly. Using which we can calculate $x$ and $y$ to get all three coordinates. The arc cosine function takes in a uniform sampling in $[-1,1]$ to give the right spread for theta

In [9]:
N=500
phi   = uniform(0,2*np.pi,size=N)
theta = np.arccos(uniform(-1,1,size=N))

st=np.sin(theta)
ct=np.cos(theta)
sp=np.sin(phi)
cp=np.cos(phi)

pts=np.c_[st*cp,st*sp,ct]

fig,ax=plt.subplots(figsize=(10,10))
ax = plt.axes(projection='3d')
ax.plot(pts[:,0],pts[:,1],pts[:,2],'o')
ax.set_aspect('equal')

# For Visualization, make spherical surface
u = np.linspace(0, 2 * np.pi, 100)
v = np.linspace(0, np.pi, 100)
x =  np.outer(np.cos(u), np.sin(v))
y =  np.outer(np.sin(u), np.sin(v))
z =  np.outer(np.ones(np.size(u)), np.cos(v))

# Plot the surface
collection=ax.plot_surface(x, y, z, color='r',alpha=0.5)


## From normal distribution¶

In this method we first get points(in 3D space) by independently drawing the x,y,z coordinates for each point from a normal distribution. The resulting points will be spherically symmetric. We then normalize the position vectors so that the points lie on a sphere. The resulting distribution will be uniform on a sphere.

In [10]:
pts=standard_normal(size=[N,3])          #get 3 N-length values
pts/=np.sqrt(np.sum(pts**2,axis=-1)[:,np.newaxis])   #normalize: v=v/|v|

fig,ax=plt.subplots(figsize=(10,10))
ax = plt.axes(projection='3d')
ax.plot(pts[:,0],pts[:,1],pts[:,2],'o')
ax.set_aspect('equal')

# For Visualization, make spherical surface
u = np.linspace(0, 2 * np.pi, 100)
v = np.linspace(0, np.pi, 100)
x =  np.outer(np.cos(u), np.sin(v))
y =  np.outer(np.sin(u), np.sin(v))
z =  np.outer(np.ones(np.size(u)), np.cos(v))

# Plot the surface
collection=ax.plot_surface(x, y, z, color='r',alpha=0.5)