= Solution
Under independent area counts and independent area <prior distributions>, the joint <posterior distribution> factors into the gamma posteriors. For each simulation $s=1,\ldots,S$, draw one joint vector by independently generating
$$
\lambda_i^{(s)}\sim\operatorname{Gamma}(y_i+1/2,E_i),\qquad i=1,\ldots,I.
$$
Rank the risks within that joint draw, taking rank one to mean highest risk:
$$
r_i^{(s)}=1+\sum_{k\ne i}\mathbf1\{\lambda_k^{(s)}>\lambda_i^{(s)}\}.
$$
The empirical <posterior rank distribution> and highest-risk probability are
$$
\boxed{\widehat{\mathbb P}(r_i=k\mid y)=\frac1S\sum_{s=1}^S\mathbf1\{r_i^{(s)}=k\},\qquad
\widehat{\mathbb P}(i\text{ highest}\mid y)=\frac1S\sum_{s=1}^S\mathbf1\{r_i^{(s)}=1\}.}
$$
Independent gamma draws make this direct <Monte Carlo integration>, so <Markov chain Monte Carlo> is unnecessary here. Continuous posteriors give ties probability zero, and the highest-risk probabilities sum to one. The same ranking calculation can be applied to dependent joint draws from a richer <hierarchical Bayesian model>; ranking separate posterior means would discard the uncertainty that the program is intended to measure.
Back to article page
