How to sample uniformly from an intersection of simplices?

Question

The title says it all. Specifically: given strictly positive real numbers $a_1,\dots,a_T$ and $b_1,\dots,b_T$, I want to sample from $$\mu:=\text{Uniform}(\{p\in[0,\infty)^T : \sum_{t=1}^T a_t p_t = 1\}\cap\{p\in[0,\infty)^T : \sum_{i=1}^T b_t p_t = 1\}),$$ assuming that the intersection is nonempty.

My best attempt so far is the following rejection sampling algorithm:

1. Set $u_t = \min\{\frac{1}{a_t},\frac{1}{b_t}\}$ for $t=1,\dots,T-1$.
2. Sample $p_t \sim \text{Uniform}(0,u_t)$ for $t=1,\dots,T-1$.
   • (Justification: if $p_t>u_t$ then $\max\{a_t,b_t\} p_t > \max\{a_t,b_t\}\min\{\frac{1}{a_t},\frac{1}{b_t}\} = 1$, so $p_t$ cannot be the $t$th coordinate of a sample from $\mu$.)
3. Set $s_a = \sum_{t=1}^{T-1} a_t p_t$ and $s_b = \sum_{t=1}^{T-1} b_t p_t$.
4. If $s_a>1$ or $s_b>1$, reject.
5. If $\frac{1}{a_T}(1-s_a) \neq \frac{1}{b_T}(1-s_b)$, reject.
6. Else, set $p_T = \frac{1}{a_T}(1-s_a)$ and accept $p=(p_1,\dots,p_T)$.

The trouble is that in step 5, we will reject with probability 1.

Does anyone know how to do this properly, either using something like the above or via a totally different approach? Thanks in advance for any ideas!

Answer 1

I think the problem lies in step 1 & 2. Setting $u_1$ in this way is fine, but once you have already sampled $p_1$, this places further constraints on what $p_2,\dots,p_T$ can be. So you should really change step 1 to be (defining $\sum_{r=1}^{0}y_r=0$):

$$u_t=\min\left(\frac{1-\sum_{r=1}^{t-1}a_rp_r}{a_t},\frac{1-\sum_{r=1}^{t-1}b_rp_r}{b_t}\right)$$

And we have a modified version of your justification at step 2:

$$max(a_t,b_t)p_t>max(a_t,b_t)u_t=\min\left(1-\sum_{r=1}^{t-1}a_rp_r,1-\sum_{r=1}^{t-1}b_rp_r\right)$$

Which leads to the same conclusion (can't be inside the sample).

And further, you have 2 constraints, so that you should not sample $p_{T-1}$ and $p_T$. So we shall "redefine" $s_a=\sum_{t=1}^{T-2}a_tp_t$ and $s_b=\sum_{t=1}^{T-2}b_tp_t$, and we have two equations:

$$1-s_a=a_{T-1}p_{T-1}+a_Tp_T$$ $$1-s_b=b_{T-1}p_{T-1}+b_Tp_T$$

Or, in matrix form:

$$ \begin{pmatrix}1-s_a\\1-s_b\end{pmatrix} = \begin{pmatrix} a_{T-1} & a_{T} \\ b_{T-1} & b_{T} \end{pmatrix} \begin{pmatrix}p_{T-1}\\p_T\end{pmatrix} $$

Which has a unique solution if $a_{T-1}b_{T}-b_{T-1}a_{T}\neq 0$. You can choose your categories for $T,T-1$ so that this condition is satisfied, because of the permutational symmetry in the problem. If you can't, then I think (but not sure) this means you really only have one constraint, and one of them can be dropped (as all of the b's are essentially functions of the a's in a similar manner to an autoregressive equation). The solution, if it exists, is given by:

$$ \begin{pmatrix}p_{T-1}\\p_T\end{pmatrix} = \begin{pmatrix} a_{T-1} & a_{T} \\ b_{T-1} & b_{T} \end{pmatrix}^{-1} \begin{pmatrix}1-s_a\\1-s_b\end{pmatrix}$$ $$=\frac{1}{a_{T-1}b_{T}-b_{T-1}a_{T}} \begin{pmatrix} b_{T} & -a_{T} \\ -b_{T-1} & a_{T-1} \end{pmatrix} \begin{pmatrix}1-s_a\\1-s_b\end{pmatrix} $$ $$\implies \begin{matrix}p_{T-1} & = & \frac{(1-s_a)b_T-(1-s_b)a_{T}}{a_{T-1}b_{T}-b_{T-1}a_{T}}\\ p_{T} & = & \frac{(1-s_a)b_{T-1}-(1-s_b)a_{T-1}}{a_{T-1}b_{T}-b_{T-1}a_{T}}\end{matrix} $$

This shows that you don't need to do rejection sampling at all for this problem