1

I'm trying to get a sample of points distributed uniformly in the 26-dimensional simplex (otherwise put, uniformly distributed 27-dimensional probability vectors). As I have been told, for example, here, the right way to do this is to sample from the uniform Dirichlet distribution.

The issue I am finding is that a largish sample from this distribution contains no points near the vertices of the simplex. In numpy:

from numpy.random import dirichlet
points = dirichlet([1] * 27, size=1000000)

After doing this repeatedly, numpy.max(points) is typically around .5. In that million samples nothing gets closer to the vertex.

I assume this is correct (and my 27-dimensional intuitions are wildly unreliable anyway). My question is: how much more should I sample in order to have, say a .9 probability of getting a point with a value of a coordinate of .9 or higher? If this turns out to be an unworkable sample, is there any other way I could coax the distribution to give me more extreme points?

Schiphol
  • 149

1 Answers1

2

Say you have $n=27$ coordinates. The probability that one of them is greater than $q=0.9$ is $n$ times the probability that a particular one of them is greater than $q$, as long as $q\ge\frac12$, since in this case these $n$ events are mutually exclusive.

Now the probability density for one coordinate $x_1$ is proportional to $\left(1-x_1\right)^{n-2}$ (since $1-x_1$ is left for the remaining $n-1$ coordinates to sum to in an $(n-2)$-dimensional simplex). The probability for $x_1\gt q$ is thus

\begin{eqnarray*} \mathsf P\left(x_1\gt q\right)&=&\frac{\int_q^1\left(1-x_1\right)^{n-2}\;\mathrm dx_1}{\int_0^1\left(1-x_1\right)^{n-2}\;\mathrm dx_1} \\ &=& \frac{\int_0^{1-q}x_1^{n-2}\;\mathrm dx_1}{\int_0^1x_1^{n-2}\;\mathrm dx_1} \\ &=& \left(1-q\right)^{n-1}\;. \end{eqnarray*}

In your case, with $n=27$ and $q=0.9$, this is $10^{-26}$, so the overall probability is $27\cdot10^{-26}$. To have a $.9$ probability for this to happen, you'd need to draw $k$ samples such that

$$ \left(1-27\cdot10^{-26}\right)^k\le1-0.9\;, $$

and solving for $k$ yields

\begin{eqnarray*} k&\ge&\frac{\log(1-0.9)}{\log\left(1-27\cdot10^{-26}\right)} \\&\simeq&-\frac1{27}10^{26}\log0.1 \\&\simeq&8.5\cdot10^{24}\;. \end{eqnarray*}

So you'd better start drawing samples right away; the universe might end any time soon.

joriki
  • 238,052
  • Thanks a lot. That's very helpful---and disheartening. Is there any way I could coax the Dirichlet distribution, or a different one, to give me more extreme points? Perhaps a way of "centering" the sample in the vicinity of the vertices? – Schiphol May 29 '18 at 11:50
  • There certainly are -- but that would then give you a different distribution. There are many ways to weight the extreme points more strongly -- it's hard to suggest one without knowing what you want to use this for. – joriki May 29 '18 at 11:52
  • 1
    I will open a new question for that, then. Thanks again. – Schiphol May 29 '18 at 11:56