There is also an analytic expression you can use, to map uniformly distributed points inside a $d$ dimensional unit hypercube, to uniformly distributed points inside a simplex of the same dimension with nodes $[\boldsymbol\xi_0, \boldsymbol\xi_1,\cdots,\boldsymbol\xi_{d+1}]$, where $\boldsymbol\xi_i\in\mathbb{R}^{d}$:
\begin{align}
M_{d}(r_1, \cdots,r_d) = \boldsymbol\xi_0 + \sum_{i=1}^{d}\prod_{j=1}^i
r_{d-j+1}^{\frac{1}{{d-j+1}}}(\boldsymbol\xi_i - \boldsymbol\xi_{i-1}),
\end{align}
The $r_{i}$ are the points inside the unit hypercube and are distributed as $\mathcal{U}[0,1]$, and the points inside the simplex (i.e. $M_d\in\mathbb{R}^d$), are also distributed uniformly. This formula is derived in an appendix of this article: https://doi.org/10.1016/j.jcp.2015.12.034
When I derived this formula I wasn't thinking of sampling on the (unit) simplex, but with a bit of experimenting, I think this is still possible by simply setting $r_d=1$, and of course by setting the nodes $\boldsymbol\xi_i$ to those of the unit simplex.
Below you'll find a Python implementation for sampling inside or on the (unit) simplex with a 2D, 3D, 4D and 100D example:
import numpy as np
import matplotlib.pyplot as plt
def sample_simplex(xi, n_mc=1, on_simplex=False):
"""
Use an analytical function map to points in the d-dimensional unit hypercube to a
d-dimensional simplex with nodes xi.
Parameters
----------
xi: array of floats, shape (d + 1, d)
The nodes of the d dimensional simplex.
n_mc : int, The default is 1.
Number of samples to draw from inside the simplex
on_simplex: boolean, default is False
If True, sample on the simplex rather than inside it
Returns: array, shape (n_mc, d)
-------
n_mc uniformly distributed points inside the d-dimensional simplex with edges xi.
"""
d = xi.shape[1]
P = np.zeros([n_mc, d])
for k in range(n_mc):
# random points inside the unit hypercube
r = np.random.rand(d)
# sample on, instead of inside the simplex (all samples will sum to one)
if on_simplex:
r[-1] = 1
# the term of the map is \xi_k_j0
sample = np.copy(xi[0])
for i in range(1, d + 1):
prod_r = 1.
# compute the product of r-terms: prod(r_{d-j+1}^{1/(d-j+1)})
for j in range(1, i + 1):
prod_r *= r[d - j]**(1. / (d - j + 1))
# compute the ith term of the sum: prod_r*(\xi_i-\xi_{i-1})
sample += prod_r * (xi[i] - xi[i - 1])
P[k, :] = sample
return P
plt.close('all')
# triangle points
xi = np.array([[0.0, 0.0], [1.0, 0.0], [0, 1]])
samples = sample_simplex(xi, 1000, on_simplex=True)
fig = plt.figure()
ax = fig.add_subplot(111)
# plot random samples
ax.plot(samples[:,0], samples[:, 1], '.')
# plot edges
ax.plot(xi[:, 0], xi[:, 1], 'o')
plt.tight_layout()
# 3D simplex
xi = np.array([[0.0, 0.0, 0.0], [1.0, 0.0, 0.0], [0, 1, 0], [0,0,1]])
samples = sample_simplex(xi, 1000, on_simplex=True)
fig = plt.figure()
ax = fig.add_subplot(111, projection='3d')
# plot random samples
ax.scatter(samples[:,0], samples[:, 1], samples[:,2], '.')
# plot edges
ax.scatter(xi[:, 0], xi[:, 1], xi[:, 2], 'o')
plt.tight_layout()
# 4D simplex
xi = np.array([[0,0,0,0], [1.0, 0.0, 0.0, 0], [0, 1, 0, 0], [0,0,1,0], [0,0,0,1]])
samples = sample_simplex(xi, 1000, on_simplex=True)
# These should all sum to one
print(np.sum(samples, axis=1))
# 100D simplex
xi = np.zeros([1, 100])
xi = np.append(xi, np.eye(100), axis=0)
samples = sample_simplex(xi, 1000, on_simplex=True)
# These should all sum to one
print(np.sum(samples, axis=1))
plt.show()
