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$

up vote 1 down vote accepted

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 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 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 at 21:18
  • Thanks @D.W. for the example. I think I know a solution... – pedroth Nov 8 at 22:51
  • 1
    @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 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 at 21:50
  • Probably this is kind of slow in high dimension (with $n = 100$). – R zu Nov 8 at 21:57

Your Answer

 
discard

By clicking "Post Your Answer", you acknowledge that you have read our updated terms of service, privacy policy and cookie policy, and that your continued use of the website is subject to these policies.

Not the answer you're looking for? Browse other questions tagged or ask your own question.