Consider identically distributed random variables $X_1, \dots, X_n$ and a constant $k \in [0,n]$. How would I sample from the distribution

$$ f_{X_1,\dots,X_n}(x_1, \dots, x_n) = c \cdot I[(\sum_{i=1}^n x_i) == k] \cdot \prod_{i=1}^n I[x_i \in [0,1]]$$

where $I$ is the indicator function and $c$ is some constant.

Essentially, I want to pick $n$ numbers in $[0,1]$ that sum to $k$ so that they are identically distributed.

Currently I am sampling $n$ realisations from a $U[0,1]$, $u_1, \dots, u_n$ and normalising to get

$$z_i = \frac{u_i \cdot k}{\sum_{i=1}^n u_i}$$

If any $z_i \notin [0,1]$ I discard the whole sample. This is currently not very efficient and I'm not sure if it is even correct (as the discarding of the normalised values may affect the distribution).

  • 1
    $\begingroup$ For special case $k=1$, a look into symmetric Dirichlet distribution may help. The marginal will be Beta distributed. $\endgroup$ – Francis Jun 17 '16 at 12:51
  • $\begingroup$ I think I have an idea. Sample a point $z$ from the symmetric Dirichlet distribution. This is essentially sampling on a simplex in $\mathbb{R}^n$ on the plane $\sum x_i = 1$. Now take one of the vertices of the simplex $s_1$ and decompose $z-s_1 = \sum_{i=2}^n a_i \cdot s_i$ to find the $a_i$. We now want point $y$ on the simplex on the plane $\sum x_i = k$. For each vertex of this simplex $r_i$ we can map to $s_i$. Therefore take $y - r_1 = \sum_{i=2}^n a_i \cdot r_i$. $\endgroup$ – rwolst Jun 18 '16 at 21:22
  • $\begingroup$ Unfortunately this doesn't look like it will work in most cases as the intersection between the hyperplane $\sum x_i = k$ and the unit hypercube is not guaranteed to be a simplex for all values of $k$. See math.stackexchange.com/questions/584006/…. $\endgroup$ – rwolst Jun 19 '16 at 7:29
  • $\begingroup$ So the intersection between the hyperplane and hypercube is a (n-1 dim) hyper simplex. Another idea would be if we knew the volume $V$ of the hypersimplex, sample $n-1$ points $u_i$ uniformly distributed on $[0,V]$. Then solve along each transformed dimension the value $r_i$ such that the section of hypersimplex with $x_i \in [0, r_i]$ has volume $u_i$. Intuitively this tells me we will have a point uniformly from the hyper simplex. However the method may not be easy in practice. $\endgroup$ – rwolst Jun 19 '16 at 7:49

A correct but inefficient method is as follows:

1) Sample $\boldsymbol{x} \in \mathbb{R}^n$ from a symmetric Dirichlet distribution (as mentioned by Francis).

2) Calculate the point $\boldsymbol{y} = k \cdot \boldsymbol{x}$.

3) If $0 \leq y_i \leq 1$ then take $\boldsymbol{y}$ as a sample else reject. The accepted samples will be distributed according to the given pdf.


As $n \rightarrow \infty$, for the worst case scenario $k$, the probability of accepting a sample generated by the algorithm $\rightarrow 0$. The argument can be made as follows:

1) The marginal probabilty of a symmetric Dirichlet distribution is $Beta(1,n-1)$.

2) In order for our sample to be accepted, we require $0 \leq x_i \leq 1/k$ for all $i = 1 \dots n$ (such that $0 \leq y_i \leq 1$).

3) We can calculate the marginal cdf to be $Pr(x_i \leq 1/k) = 1 - (1 - 1/k)^{n-1}$.

4) Choose $k = n/c$ for some integer $c \geq 2$. Note we do not consider values of $k$ greater than $n/2$ as by symmetry they are equivalent to a value less than $n/2$.

5) $Pr(x_i \leq c/n) = 1 - (1-c/n)^{n-1} \rightarrow 1-e^{-c}$ as $n \rightarrow \infty$.

6) Therefore take $c=2$ for the worst case scenario.

7) (Hand waving) As $n \rightarrow \infty$, for $k = n/2$, for some $a n^b < n_0 < n$ where $a>0$ and $0 < b < 1$, the joint distributions of $n_0$ of our $x_i$ become almost independent.

8) Therefore $Pr(x_1, \dots, n_{n_0} \leq 1/k) \approx \prod_{i=1}^{n_0}Pr(x_i \leq1/k) = (1-e^{-2})^{n_0}$

9) So the probability of $n_0$ of the samples being less than $1/k$ tends to $0$ as $n_0 \rightarrow \infty$.

10) Therefore the probability the method chooses a sample $\boldsymbol{x}$ such that all $x_i \leq 1/k$ tends to 0 as $n \rightarrow \infty$.


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.