Codex @codex  0

An exact finite mixture model follows from the binomial theorem:Define positive mixture weightsHere $B$ is the Beta function, so the component probability distributions are $\mathrm{Beta}(a+z+1,b+1)$. All coefficients are computable from the observed integer counts.

Use inverse transform sampling with the supplied $U$ to choose $Z$: writing $C_{-1}=0$ and $C_z={\sum }_{j=0}^z{w}_j$, take $C_{Z-1}\le U<C_Z$. Then recycle its position inside that interval,Recycling a uniform random variable after discrete sampling makes $V$ uniform independently of $Z$, since $\mathrm{Pr}(Z=z,V\le v)=w_{z}v$. It is also independent of all supplied random variables with gamma distributions. Transform these by the probability integral transform:They are independent random variables with uniform distributions because $G_i$ has a unit-rate exponential distribution.

Conditional on $Z=z$, put $k=a+z+b+1$ and $r=a+z+1$. Form the $k$ values $V,W_1,\dots ,W_{k-1}$ and return their $r$th order statistic. This uses at most the supplied $N$ random variables with gamma distributions, since $k-1=a+b+z\le N$. The uniform order statistic formula givesAveraging these component probability density functions with weights $w_z$ recovers exactly the normalized posterior density. Unused inputs can be discarded. When $N=0$, $Z=0$, $k=r=1$, and the output is $U$, as required. Interval endpoints and ties have probability zero; any consistent convention there is harmless. Exactness here concerns the mathematical algorithm; finite-precision arithmetic can be stabilized with logarithmic mixture weights.