# Uniform sampling over non-standard simplex

Uniform sampling over a $$n$$-dimensional standard simplex is described here: Uniform sampling from a simplex

I want to sample one point from a non-standard simplex with vertices at:

• $$s_{i}\vec{e_{i}}$$, where $$\vec{e_{i}}$$ is the $$i^{th}$$ vector in the standard basis, and $$s_{i} \ge 0$$.
• the origin is : $$\vec{0}$$

How to do that? It seems to be no longer a special case of the Dirichlet distribution.

Hit-and-run sampling would need a lot of cycle to diffuse throughout the simplex. When I only need one point, the other points are wasted.

Currently, $$n \approx 15$$

According to Wikipedia the $$n$$-simplex is the set defined by:

$$S = \left \{ \mathbf{x} \in \mathbb{R}^n \; | \; \mathbf{x} = \sum_{i=0}^n \mathbf{u}_i\theta_i , \; \mathbf{u}_k \in \mathbb{R}^n , \;\; \theta_k > 0 \;\forall k, \;\;\sum_{i=0}^n \theta_i = 1\right \}$$

Note that the $$\boldsymbol{\theta} = (\theta_0, \dots, \theta_n)^T$$ is the standard $$n$$-simplex.

Hence in order to sample from $$S$$, compute a linear combination of the vertices $$\mathbf{u}_k$$ with the weights given by a sample of the standard simplex.

Which is the same as the following random process:

Consider the matrix of vertices $$U = [\mathbf{u}_0 \; \dots \; \mathbf{u}_n]$$, then a sample in $$S$$ is given by:

$$\boldsymbol{\theta} \sim \text{Standard-Simplex}$$ $$\mathbf{x} = U \boldsymbol{\theta}$$

Note that in your case the vertex $$\mathbf{u}_k = s_k \mathbf{e}_k$$.

• In order to be a scaled version, all the elements of the diagonal should be equal, which in your case is not. – pedroth Nov 8 '18 at 17:21
• I really really hope this works. Please tell me if it does. The argument that this doesn't work is: If scaling works, I can scale several axes to zero (set several $s_k = 0$). But then, when I scale a triangle to a line, the distribution on the line is no longer uniform. Therefore, it doesn't work. – R zu Nov 8 '18 at 17:28
• This answer is incorrect. It doesn't lead to a uniform distribution on $\mathbf{x}$. As a simple example, consider $U = \begin{pmatrix} 2 &0\\ 0&1 \end{pmatrix}$; then if $\boldsymbol{\theta}$ is uniformly distributed (on some set), $\mathbf{x} = U\boldsymbol{\theta}$ won't be uniformly distributed (on a corresponding set). – D.W. Nov 8 '18 at 21:18
• Thanks @D.W. for the example. I think I know a solution... – pedroth Nov 8 '18 at 22:51
• @pedroth Sorry. I think you are right. This actually works if all $s_{i} \ne 0$ becuase $U$ is a linear transformation that is not singular. That means the determinant of $U$ is constant and nonzero. By change of coordinate formula in calculus, the probability density on the non-standard simplex is scaled by the determinant of $U$, which is same for all the points on the standard simplex. – R zu Nov 9 '18 at 16:30

I make this up for $$n$$-dimensional simplex.

Looks ok but I am not sure.

Barycentric random walk

1. Choose an initial point $$\vec{x}$$
2. Choose $$\lambda \sim (U[0, 1])^{1/n}$$
3. Choose a random vertex of the simplex $$\vec{v}$$
4. $$\vec{x} \leftarrow \lambda\vec{x} + (1 - \lambda)\vec{v}$$
5. Repeat till $$\vec{x}$$ is sufficiently random.

Justification

The probability density of $$\lambda$$ is proportional to the cross sections along the line between the current point and the vertex.

Code

import numpy as np

N = 10000

scales = np.array([0, 200, 3])

# -- Barycentric random walk --
vertices = np.diag(scales)
nd = scales.size
x = scales.copy().astype(float)

x[1:] = 0

x_all = np.empty((N, nd))
for i in range(N):
# choose random vertex
t = np.random.uniform(0, 1.0) ** (1 / nd)
j = np.random.choice(nd)
x = x * t + vertices[j] * (1 - t)
x_all[i] = x

x = x_all

# -- Plot --

from mpl_toolkits.mplot3d import Axes3D
import matplotlib.pyplot as plt

fig = plt.figure()
ax = fig.add_subplot(111, projection='3d')

xs, ys, zs = x[:, 0], x[:, 1], x[:, 2]
ax.scatter(xs, ys, zs, marker='.')
ax.set_xlabel('X')
ax.set_ylabel('Y')
ax.set_zlabel('Z')

plt.show()

• Are there any random walk algorithm over a simplex in the literature? – R zu Nov 8 '18 at 21:50
• Probably this is kind of slow in high dimension (with $n = 100$). – R zu Nov 8 '18 at 21:57