# 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. Jun 17, 2016 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$. Jun 18, 2016 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/…. Jun 19, 2016 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. Jun 19, 2016 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$.

You can use the hitandrun R package to generate points uniformly over some regions. As an example from your question, let's say you want $$X_1 + X_2+X_3+X_4=2,$$ where $$0 \leq X_1, X_2, X_3,X_4 \leq 1$$

Here is the code to do this, where I've included checks so you can assess whether indeed the components are identically distributed:

library(hitandrun)
library(data.table)

n <- 4

# Here we force the sum of the 4 variables to be 2
eq.constr <- list(
constr = matrix(c(1,1,1,1),nrow=1,ncol=n),
dir = '=',
rhs = 2)

# The following bounds the variables between 0 and 1
ineq.constr <- list(
constr = rbind(-diag(n),diag(n)),
dir = rep('<=',n*2),
rhs = c(rep(0,n),rep(1,n)))

basis <- solution.basis(eq.constr)
transform <- createTransform(basis)
constr <- transformConstraints(transform,ineq.constr)
x0 <- createSeedPoint(constr,homogeneous=TRUE)
x <- har(x0,constr,100000,transform=transform, homogeneous=TRUE)\$samples

stopifnot(all.equal(x[,1],x[,2],x[,3],2-x[,4]))
stopifnot(all(x >= 0))
stopifnot(all(x <= 1))

# Thin the output to lessen correlation
x <- x[seq(10,nrow(x),10),]

# Are the component variables identically distributed?
ks.test(x[,1],x[,2])
ks.test(x[,1],x[,3])
ks.test(x[,1],x[,4])
ks.test(x[,2],x[,3])
ks.test(x[,2],x[,4])
ks.test(x[,3],x[,4])

# Get a histogram of a marginal distribution
hist(x[,1])


Notes: I did have some issues with large initial samples (for example, a million variates). So you'll have to experiment to see what works for you. Similarly, I chose the thinning of using every 10th observation arbitrarily - you'll want to do your own analysis of what is needed for your situation.

Here is some sample output of the variates:

           V1        V2        V3         V4 total
1: 0.25230984 0.5435753 0.8445336 0.35958126     2
2: 0.08133462 0.7457540 0.9077554 0.26515597     2
3: 0.10021822 0.4472125 0.5515221 0.90104719     2
4: 0.80093829 0.3557503 0.8008299 0.04248152     2
5: 0.09021195 0.3081201 0.7329975 0.86867048     2
6: 0.48967179 0.1687135 0.9135057 0.42810901     2