# Generating identically-distributed random variables with a constraint

Is there a way to generate identically-distributed random variables (eg $x_1,x_2,x_3,x_4$) with the following constraint:

$\frac{x_1*x_2}{x_3*x_4} ≡ 1$

$x \in (0,1)$

Please note that simply sampling $x_1, x_2, x_3$ from the given distribution and then finding $x_4$ as $\frac{x_1*x_2}{x_3}$ would not work since $x_4$, in that case, would not follow the given distribution.

Edit1: One inefficient way of doing it may be generating a large number of $x_1,x_2,x_3,x_4$ and then simply filtering with some tolerance the ones satisfying the constraint. The method may be not as inefficient here because of this specific constraint (ratio = 1) since we are going to get lots of instances where the constraint is approximately satisfied.

• Are these variables necessarily positive? Is anything else known/specified? Is any kind of dependence permissible, or are there limits on the amount and kind of dependence? – Glen_b -Reinstate Monica May 13 '18 at 8:01
• without having other information, I suggest to find the distribution of $z=(x_1x_2)/(x_3x_4)$,under the constraint, that is easy to find, and to use one of the many sampling algorithm – niandra82 May 13 '18 at 10:13
• Glen_b, yes, there are some specifics. x is a positive number between 0 and 1. I will add it in the post. As of now I try not to limit the variable to any particular distribution. – PG_eon May 14 '18 at 0:45
• Niandra, thank you for your suggestion. Could you elaborate a bit on what you mean by "finding ... under constraint"? I am not sure how what you are suggesting is different from my Edit1 – PG_eon May 14 '18 at 0:48
• Your constraint is $x_1x_2=x_3x_4$. So first draw $x_1$ and $x_2$ independently. Then form their product $x_1x_2$. Then draw $x_3$ independently. If $d = x_1x_2/x_3 < 1$, then set $x_4 = d$ and you have completed one sample. Otherwise, reject $x_3$ and independently re-sample $x_3$ until it satisfies the criterion, and then set $x_4 = d$, and you have completed the sample. I leave it to you as an exercise to determine whether this method generates samples from the correct conditional distribution, or whether I have screwed up. – Mark L. Stone May 14 '18 at 1:48

In your problem you wish to generate a set of four random variables $X_1, X_2, X_3, X_4$ with identical marginal distributions that have support on the interval $0<x<1$, and subject to the additional constraint $X_1 \cdot X_2 = X_3 \cdot X_4$. We will refer to the marginal distribution of each of these values as the desired marginal. You can generate random variables obeying these constraints as follows.

Preliminary: Start by choosing any bivariate symmetric joint density function $p_{A,B}$ that has the desired marginal densities for both variables. From this joint density, you can define $D = A \cdot B$ and you obtain the corresponding joint density:

$$p_{A, D}(a, d) = \begin{vmatrix} 1 & 0 \\ d/a & a \end{vmatrix} p_{A,B}(a, d/a) = a \cdot p_{A,B}(a, d/a) \quad \quad \text{for all }0<d<a.$$

You then have the corresponding densities:

$$p_{A| D}(a, d) = \frac{p_{A,D}(a, d)}{p_{D}(d)} \quad \quad \quad p_{D}(d) = \int_d^1 a \cdot p_{A,B}(a, d/a) da.$$

Generating your values: Now, once you have got these density functions, you are ready to generate your sample as follows:

• Generate a single value $D \sim p_D$;
• Now use this value to generate $X_1, X_3 | D \sim \text{IID }p_{A|D}$;
• Set $X_2 = D / X_1$ and $X_4 = D / X_3$.

This will give you generated random variables $X_1, X_2, X_3, X_4$ with marginal distribution $p_A$ (which is the desired distribution) and subject to the constraint $X_1 \cdot X_2 = X_3 \cdot X_4$. Note that the method works for any choice of symmetric bivariate distribution $p_{A,B}$ and so the solution to the problem is non-unique. The particular choice of $p_{A,B}$ determines the dependency structure.

Special case: A simple special case is to choose a starting bivariate joint density where the values are independent with the specified desired marginal distribution. This gives you the simplified form $p_{A,B}(a,b) = p_A(a) \cdot p_A(b)$, and the above equations simplify to:

$$p_{A, D}(a, d) = a \cdot p_{A}(a) p_A(d/a) \quad \quad \text{for all }0<d<a,$$

and:

$$p_{A| D}(a, d) = \frac{a \cdot p_{A}(a) p_A(d/a)}{p_{D}(d)} \quad \quad \quad p_{D}(d) = \int_d^1 a \cdot p_{A}(a) p_A(d/a) da.$$

Given a desired distribution $p_A$, it should not be difficult to derive these distributions and perform the generation algorithm. (Note that even though you start with a bivariate density with independence, you still obtain generated values that are not independent, since they obey your linear constraint.)

• Your scheme is not really very different from what I outlined in the comments on the original post, except I formed the joint distribution by simulation rather than analytically. However, you have neglected to account for the possibility of the need to reject some samples, which you would be subject to, because $D/X_1$ or $D/X_3$ may be out of bounds. – Mark L. Stone May 14 '18 at 2:29
• The analytic method specified always gives $D < X_1 < 1$, since this is the support of the conditional distribution. Hence, the value $D/X_1$ should always be within the required bounds (i.e., there is no need for any rejection step). (The same applies for $X_3$ and $D/X_3$.) – Ben - Reinstate Monica May 14 '18 at 2:51
• @Ben Perhaps you might know the answer to the following similar question? stats.stackexchange.com/questions/390049/… – Kagaratsch Jan 30 '19 at 23:05

Let's look at a special case and use natural logs to convert this to a sum. Suppose we have $X_1,X_2 \sim U[0,1]$ with $X_1 X_2 = z.$ Let $R= - \mathrm{ln} \left( X_3 \right)$ and $S= - \mathrm{ln} \left( X_4 \right).$

Then we require $R + S = - \mathrm{ln} \left( X_3 \right) + - \mathrm{ln} \left( X_4 \right) = - \mathrm{ln} \left( z \right).$

It can be determined that $R$ and $S$ are gamma random variables with both parameters equal to 1 (exponential in this case).

So now we have a fixed sum of iid gamma random variables. Applying the general logic outlined in the answer here to the gamma case (How to generate two groups of $n$ random numbers in $U(0,1)$ such that sum of these two groups equal?), we find the conditional distribution of $R$ given $z$ is $$R|z \ \sim \ U \left[ 0,-\mathrm{ln} \left( z \right ) \right]$$

Then the algorithm is:

1. Generate independent realizations $x_1$ and $x_2$ as $U[0,1].$
2. Calculate $z=x_1x_2.$
3. Generate a realization $r$ from $R \sim U[0,- \mathrm{ln} \left( z \right) ]$
4. Calculate $x_3 = e^{-r}$
5. Calculate $x_4 =z/x_3$

On taking a second look, we can write $-r = \left[ {\mathrm{ln}} \left( z \right) \right] U,$ where $U$ is uniform on $[0,1].$ Then $e^{-r} = z^U.$

Combining some of the above steps, we can write the process as

1. Generate $x_1,x_2$ as independent $U[0,1]$ realizations and find $z=x_1x_2.$
2. Generate $U \sim U[0,1]$ and calculate $x_3 = z^U$ and $x_4=z^{\left( 1 - U \right)}=z/x_3$

This is an extended comment asking for clarification of @MarkL.Stone 's comment. If I following the procedure correctly, I don't get x1, x2, x3, and x4 having the same marginal distributions if I start out with x1, x2, and x3 having uniform distributions on (0,1):

n <- 10000
x1 <- runif(n)
x2 <- runif(n)
x3 <- rep(NA,n)
for (i in 1:n) {
x3[i] = runif(1)
while ((x1[i]*x2[i]) >= x3[i]) {
x3[i] <- runif(1)
}
}
x4 <- x1*x2/x3
par(mfrow=c(2,2))
hist(x1, freq=FALSE, main="x1", ylim=c(0,2.5), las=1)
hist(x2, freq=FALSE, main="x2", ylim=c(0,2.5), las=1)
hist(x3, freq=FALSE, main="x3", ylim=c(0,2.5), las=1)
hist(x4, freq=FALSE, main="x4", ylim=c(0,2.5), las=1)


What did I interpret incorrectly?