 Solution

In Markov chain Monte Carlo, retain the derived quantity $\delta =p_M-p_T$ at every iteration. Rough BUGS code using the actual observations isSupply milkAnswers=3 and teaAnswers=1. Summarize delta by its posterior mean, empirical quantiles and credible interval; the average of positive estimates $\mathbb{P}(p_M>p_T\mid \mathcal{D})$. Here $\mathbb{E}[\delta \mid \mathcal{D}]=0.4$. Check Markov chain Monte Carlo convergence diagnostics before interpreting the simulation. Since both Bayesian posteriors are independent known Beta distributions, direct independent sampling is an equally valid, simpler way to obtain the same summaries.