Given a fixed positive integer $n$, I have two random variables $$A(n)=2^{A\_2}\cdots p^{A\_{p\_n}}, B(n)=2^{B\_2}\cdots p^{B\_{p\_n}},$$ where $p\_n$ is the largest prime number not exceeding $n$, $(A\_p)\_{p\le n}$ and $(B\_p)\_{p\le n}$ are two random processes, and the index $p$ ranges over prime numbers. I am searching for a coupling of $A(n), B(n)$ such that $A\vert B$, and this holds if and only if we have the vector inequality $(A\_2,A\_3\ldots, A\_{p\_n})\le (B\_2,B\_3\ldots, B\_{p\_n})$.


Thus, I seek a probability space in which $A\_p\le B\_p$ for each $p\le n$. 


For each $p$, I have a coupling of $A\_p$ and $B\_p$ such that $A\_p\le B\_p$. I.e., for each $p$ I have a joint distribution $(A\_p',B\_p')$ with marginals $A\_p,B\_p$ such that $\mathbb{P}(A\_p=a, B\_p=b)=0$ when $a > b$.


Is there any literature/technique on how to combine these couplings into a new probability space in which $A\_p\le B\_p$ for *all* $p\le n$?


My idea was to consider the joint distribution of the $n$ coupled random vectors $(A\_p,B\_p)$. Then the $p$th marginal is $(A\_p, B\_p)$ which only has positive probability when $A\_p \le B\_p$. However, I am not sure if this is the desired probability space since the marginals are not the $A\_p$'s or $B\_p$'s but are instead random vectors. I am unsure if this joint distribution gives the desired coupling of $A,B.$



Combining Couplings of Random Variables
One way to achieve your goal is as follows. You have nonnegative numbers $f\_p(a\_p,b\_p)=P(A\_p=a\_p,B\_p=b\_p)$ for $p=2,3,5,\dots$, $a\_p=0,1,\dots$, $b\_p=0,1,\dots$ such that for each $p$ you have $\sum\_{a\_p,b\_p}f\_p(a\_p,b\_p)=1$ and $f\_p(a\_p,b\_p)=0$ when $a\_p>b\_p$. Consider the corresponding conditional probabilities $f\_p(a\_p|b\_p):=f\_p(a\_p,b\_p)/g\_p(b\_p)$, where $g\_p(b):=\sum\_{a\_p}f\_p(a\_p,b)=P(B\_p=b)$, your marginal probabilities for $B\_p$ -- assuming $g\_p(b)>0$ for all relevant $b$. You can use these couplings of $A\_p$ and $B\_p$ by requiring that the sequence $(A\_2,A\_3,A\_5,\dots)$ be conditionally independent given $(B\_2,B\_3,B\_5,\dots)$ and that each $A\_p$ depend on the sequence $(B\_2,B\_3,B\_5,\dots)$ only through $B\_p$: 
\begin{equation}
 P(A\_2=a\_2,A\_3=a\_3,A\_5=a\_5,\dots|B\_2=b\_2,B\_3=b\_3,B\_5=b\_5,\dots)
 =f\_2(a\_2|b\_2)f\_3(a\_3|b\_3)f\_5(a\_5|b\_5)\dots
\end{equation}
for all $a\_p$'s and $b\_p$'s. For the joint distribution of $(B\_2,B\_3,B\_5,\dots)$ you can take any distribution with the marginals $(g\_2,g\_3,g\_5,\dots)$. It is easy to check that thus you will have a joint distribution of $(A\_2,A\_3,A\_5,\dots,B\_2,B\_3,B\_5,\dots)$ with the prescribed two-dimensional marginals $f\_p$: $P(A\_p=a\_p,B\_p=b\_p)=f\_p(a\_p,b\_p)$ for all $p=2,3,5,\dots$, $a\_p=0,1,\dots$, $b\_p=0,1,\dots$. 


In particular, letting $B\_2,B\_3,B\_5,\dots$ be independent, you will just get independent pairs $(A\_p,B\_p)$ with the prescribed two-dimensional marginals $f\_p$. 


