Solution (source code)

= Solution

Put $a=y_4$ and $b=y_2+y_3$. Dropping factors independent of $\theta$, the <multinomial likelihood> and <uniform prior> give the <posterior density>
$$
\pi(\theta\mid y)\propto(2+\theta)^{y_1}\theta^a(1-\theta)^b\mathbf1_{(0,1)}(\theta).
$$
Multiplying this by the stipulated <binomial distribution> for the <latent variable> cancels the factor $(2+\theta)^{y_1}$:
$$
\pi(\theta,z\mid y)\propto\binom{y_1}{z}2^{y_1-z}\theta^{a+z}(1-\theta)^b,
\qquad 0\leq z\leq y_1.
$$
Consequently the two <full conditional distributions> are
$$
\boxed{Z\mid\theta,y\sim\operatorname{Binomial}\!\left(y_1,\frac{\theta}{2+\theta}\right),\qquad
\theta\mid Z=z,y\sim\operatorname{Beta}(a+z+1,b+1).}
$$
Start with $0<\theta_0<1$, sample $Z$ from its <binomial distribution>, then sample a new $\theta$ from its <Beta distribution>, and repeat. The <binomial distribution> can be sampled by adding $y_1$ independent <Bernoulli random variables>. For the <Beta distribution>, take independent $A\sim\operatorname{Gamma}(a+z+1,1)$ and $B\sim\operatorname{Gamma}(b+1,1)$ and return $A/(A+B)$; integer-shape <gamma distributions> are sums of independent unit-rate <exponential distributions>. Both shapes are positive even when some observed counts vanish.

Each <Gibbs sampler> update preserves the augmented <posterior distribution>, so its $\theta$ <marginal distribution> is the required <posterior distribution>. The positive <full conditional distributions> on the interior allow exploration of the entire support. These iterates are generally dependent; they are not the independent exact draws constructed in the later parts.