2 Random 2 Sampling

Abhilash Reddy M, May 28 2022

This is a sequel to my previous Jupyter notebook on random sampling where I covered generating uniformly random points on various domains. This notebook demonstrates sampling in some additional domains and some additional methods. 3D plots herein are drawn using plotly and are interactive. Some online notebook viewer websites do not support interactive plots. So, the plotly plots might not appear. This notebook is available here. For best the experience please read notebook hosted here or run this notebook somewhere that supports all the features and modules used.

Contents

1. Simplex
   • Triangles
   • Tetrahedra
   • An exercise
2. Sphere

We begin by setting up the notebook environment

%matplotlib inline
from jupyterthemes import jtplot
jtplot.style() # set the style of matplotlib plots to match the notebook theme
import numpy as np
from numpy.random import uniform,standard_normal
import matplotlib.pyplot as plt
import plotly.graph_objs as go
import plotly.io as pio
pio.templates.default = "plotly_dark" #make plotly use the plotly_dark theme

Sampling in a simplex

A Simplex is a generalization of triangles (2D) and tetrahedra (3D) to arbitrary dimensions. The name is because it refers to the simplest polytope in a given space. There is a great general algorithm that generates random points in any simplex.

Sampling in an $K$-dimensional simplex

1. Let $V_{ij}$ be the $[(K+1)\times K]$ array of the vertex coordinates that form a $K$-d simplex
2. Let $u_i$ be $K+1$ uniform random samples from (0,1]
3. Let $v_i = \dfrac{\ln{u_i}}{\sum_{j=1}^{K+1}{\ln{u_j}}}$
4. $p_i = V_{ji}v_j$ is uniformly random within the simplex

Thats it!

(I have used index notation above. In $V_{ij}$ $i$ refers to the vertex and $j$ refers to coordinate dimension)

We have started with samples from a uniform distribution and in the second step converted them to an exponential distribution (actually $exp\sim-\ln{u}$, but since we normalize, the signs cancel and hence I have left out the negative signs in the algorithm above) and normalized so that their sum is unity. Now, these are valid barycentric coordinates for some point within the simplex. The random point is simply the weighted sum of the vertices with the barycentric coordinates as weights.

This method needs one extra uniform sample compared so some other methods for lines and triangles. Within a line we can specify a unique point with just one coordinate and within a triangle two barycentric coordinates are enough to identify a unique point. This is because the sum of the barycentric coordinates should be unity for any point inside the simplex. This is not a comment on the efficiency of any method. In fact, I will not comment on the efficiency of the methods here.

A demonstration of this algorithm in 2D is shown below

triverts=np.asarray([[0,0],[1.0,0.2],[0.4,0.9]])
fig,ax=plt.subplots(figsize=[4,4], subplot_kw=dict(aspect='equal',))
_=ax.add_artist(plt.Polygon(triverts,alpha=0.3,lw=2,ec='k'))
N=1000
pts=[]
for i in range(N):
    pqr=uniform(size=3)
    pqr=np.log(pqr)
    pqr=pqr/np.sum(pqr)
    pts.append(triverts.T.dot(pqr))
pts=np.asarray(pts)
_=ax.plot(pts[:,0],pts[:,1], 'ro', alpha=0.5)

We can write a general dimension-agnostic function like so:

def simpsRUS(vertices,N):
    """
    Random Uniform Sampling in a K-d simplex
    Input:
        vertices: a list of vertices of the simplex
        N: number of points to be generated
    """
    pts=np.empty((N,vertices.shape[1]))
    n=len(vertices)
    for i in range(N):
        pqr=uniform(size=n)
        pqr=-np.log(pqr)
        pqr=pqr/np.sum(pqr)
        pts[i,:]=vertices.T.dot(pqr)
    return pts

The following block uses this function to get random points inside a triangle

triverts=np.asarray([[0,0.4],[1.1,0.3],[0.4,1.1]])
fig,ax=plt.subplots(figsize=[4,4], subplot_kw=dict(aspect='equal',))
_=ax.add_artist(plt.Polygon(triverts,alpha=0.3,lw=2,ec='k'))
N=1000
pts=simpsRUS(triverts,N)
_=ax.plot(pts[:,0],pts[:,1], 'ro', alpha=0.5)

Tetrahedron

Lets try to use this function to sample points inside a tetrahedron. I am going to use plotly to make interactive 3D plots

# we need to make a trimesh to plot the tetrahedron
tetverts=np.asarray([[0.1,0,0],[1.1,0.1,0],[0,1.3,0],[1,0,1.2]])
triangles=np.asarray([[0,1,2],[0,1,3],[0,2,3],[1,2,3]])

x, y, z = tetverts.T
i,j,k = triangles.T

N=2000
pts=simpsRUS(tetverts,N)

fig = go.Figure(data=[go.Mesh3d(x=x, y=y, z=z, i=i, j=j, k=k, color='lightblue', opacity=0.25,)])
fig.add_trace(go.Scatter3d(x=pts[:,0], y=pts[:,1], z=pts[:,2], mode='markers', marker=dict(size=2, color='red', opacity=0.5)))
fig.show()