3
$\begingroup$

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

$\endgroup$
  • 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
1
$\begingroup$

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

| cite | improve this answer | |
$\endgroup$
0
$\begingroup$

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

| cite | improve this answer | |
$\endgroup$

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.