9
$\begingroup$

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.

Improve this question
$\endgroup$
3
  • 4
    $\begingroup$ Simultaneously corssposted to cs.se - please don't do this. cs.stackexchange.com/questions/60311/… $\endgroup$
    – usul
    Jul 5, 2016 at 3:17
  • 2
    $\begingroup$ Your solution in Edit 2 is the right algorithm to use: you have identified the vertices of the simplex, now just generate a uniform point on the probability simplex (using one of the various known methods for this), and use the sampled point to assign coefficients to the vertices and output the resulting convex combination. I suggest you post the solution yourself and accept it. $\endgroup$
    – Sasho Nikolov
    Jul 14, 2016 at 0:28
  • 1
    $\begingroup$ Oh, yes, that is exactly what I meant. I did not post on the algorithm on drawing from the prob simplex as that is already solved in one of the links I posted. $\endgroup$
    – Juanito
    Jul 14, 2016 at 0:55

1 Answer 1

Reset to default
2
$\begingroup$

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

Improve this answer
$\endgroup$
3
  • 6
    $\begingroup$ These algorithms are quite the overkill for the problem. Sampling from any simplex $\Delta$ reduces to sampling from the probability simplex as long as you can identify the vertices of $\Delta$. Identifying the vertices in this case is an exercise. There are a number of (simple) classical algorithms for sampling from the probability simplex. $\endgroup$
    – Sasho Nikolov
    Jul 14, 2016 at 0:32
  • $\begingroup$ @SashoNikolov Is it ok if multiple $a_i$ are zero? In that case, the transformation from the standard simplex to non-standard simplex is linear and singular. Projection of a uniform distribution of points over an equilateral triangle to the base would give a non-uniform distribution of points on the base. The center of the base will have more points. $\endgroup$
    – R zu
    Nov 9, 2018 at 13:56
  • $\begingroup$ @Rzu Any polytope is a (possibly singular) linear transformation of the probability simplex. Sampling from an arbitrary polytope given by its vertices is more complicated, I believe, and probably needs the full power of Monte Carlo methods. But in this case, I don't see why you say that the linear transformation is singular. $\endgroup$
    – Sasho Nikolov
    Nov 10, 2018 at 23:59

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.