3
$\begingroup$

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$

$\endgroup$
0

2 Answers 2

2
$\begingroup$

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$.

$\endgroup$
5
  • $\begingroup$ In order to be a scaled version, all the elements of the diagonal should be equal, which in your case is not. $\endgroup$
    – pedroth
    Nov 8, 2018 at 17:21
  • $\begingroup$ 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. $\endgroup$
    – R zu
    Nov 8, 2018 at 17:28
  • $\begingroup$ 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). $\endgroup$
    – D.W.
    Nov 8, 2018 at 21:18
  • $\begingroup$ Thanks @D.W. for the example. I think I know a solution... $\endgroup$
    – pedroth
    Nov 8, 2018 at 22:51
  • 1
    $\begingroup$ @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. $\endgroup$
    – R zu
    Nov 9, 2018 at 16:30
0
$\begingroup$

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()
$\endgroup$
2
  • $\begingroup$ Are there any random walk algorithm over a simplex in the literature? $\endgroup$
    – R zu
    Nov 8, 2018 at 21:50
  • $\begingroup$ Probably this is kind of slow in high dimension (with $n = 100$). $\endgroup$
    – R zu
    Nov 8, 2018 at 21:57

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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