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.

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