Generating identically-distributed random variables with a constraint

Question

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.

Answer 1

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

Answer 2

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$

Answer 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?