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.
  • 167,959
  • 22
  • 232
  • 500
pedroth
  • 245
  • 1
  • 8
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