[Math] Uniform sampling from general simplex with a twist

algorithmslinear programmingmg.metric-geometrypr.probabilityprobability distributions

This is part of a question I had asked elsewhere, and then some of the links redirected me to CS stack exchange.

Given $0\leq a_1\leq\dots\leq a_D\leq1$ (all strictly positive), I want to draw points uniformly from

$\qquad X =\{(x_1, x_2, …, x_D) \mid 0 \leq x_D \leq x_{D-1} \leq \dots \leq x_1\ , \sum_{i=1}^D a_i x_i = 1 \} $

How do I achieve this? I would also like to understand how $x_1$ might be distributed over $[\frac {1} {\sum_i a_i}, \frac {1} {a_1}$] (a theoretical query).

This question solves it for $a_i=1$ and without the $x_i$'s having a prior order, but I do not see any obvious way to generalize from that answer.

There is a similar question on MO, but without the weights $a_1,\dots,a_D$.

Edit for possible solution clue:

I think this page from Devroy[1] page 568 is really helpful. To the best of my understanding, it solves my problem except for the part of $x_1 \geq x_2 \geq \dots \geq x_D$.

enter image description here

I think I might have an idea of how to solve it all the way, and I will try to write it out in the next couple of days.

Edit 2: To complete the answer for $x_1 \geq x_2 \geq \dots \geq x_D$ part, one needs to note that for any $(x_1,x_2,x_3)$,
$$(x_1,x_2,x_3)=(\frac {1} {a_1},0,0)*(x_1-x_2)a_1+(\frac {1} {a_1+a_2},\frac {1} {a_1+a_2},0)*(x_2-x_3)(a_1+a_2)+(\frac {1} {a_1+a_2+a_3},\frac {1} {a_1+a_2+a_3},\frac {1} {a_1+a_2+a_3})*x_3*(a_1+a_2+a_3)$$

In fact this representation is an if and only if, as in, any convex combination of the said points will be in the aforementioned plain $X$ and would have $x_1 \geq x_2 \geq \dots \geq x_D$.

And this can be extended to a general number of dimensions.

So to draw points uniformly from $X$ with $0 \leq x_1 \leq x_2 \leq \dots \leq x_D$, one needs to start by drawing numbers uniformly on the unit simplex, and then use that draw as weights to find a convex combination of the points
$(\frac {1} {a_1},0,0), (\frac {1} {a_1+a_2},\frac {1} {a_1+a_2},0), (\frac {1} {a_1+a_2+a_3},\frac {1} {a_1+a_2+a_3},\frac {1} {a_1+a_2+a_3})$. Similarly for higher dimensions.

Best Answer

I believe there is something called the Metropolis Algorithm The volume of your simplex (without the linear constraint) is $\frac{1}{n!}$ and so rejection sampling is terrible.

So instead you do a reflecting random walk in your space and the limiting distribution is uniform in the polyhedron.

If you want to know rigorously why the Metropolis-Hastings algorithm works, maybe you can try What do we know about the Metropolis-Hastings algorithm? by Diaconis and Salff-Lacoste.

Oh! I found it... Gibbs/Metropolis algorithms on a convex polytope ... but these very complicated. Just remember the phrase "detailed-balance".

Related Question