3

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$

R zu
  • 168
  • 10

2 Answers2

2

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

D.W.
  • 159,275
  • 20
  • 227
  • 470
pedroth
  • 235
  • 1
  • 8
  • In order to be a scaled version, all the elements of the diagonal should be equal, which in your case is not. – pedroth Nov 08 '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 08 '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 08 '18 at 21:18
  • Thanks @D.W. for the example. I think I know a solution... – pedroth Nov 08 '18 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 09 '18 at 16:30
0

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()
R zu
  • 168
  • 10