# Sample From a 'Constrained' Distribution

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

• For special case $k=1$, a look into symmetric Dirichlet distribution may help. The marginal will be Beta distributed. – Francis Jun 17 '16 at 12:51
• 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$. – rwolst Jun 18 '16 at 21:22
• 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/…. – rwolst Jun 19 '16 at 7:29
• 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. – 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.

Inefficency

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