In the case of comparing two or more populations, we might be interested in the posterior distribution of $|\theta_1-\theta_2|, \; \theta_1/\theta_2,$ or { $\;\theta_1,...,\theta_m$}, which are functions of more than one parameter.
Obtaining exact values for these posterior quantities can be difficult or impossible, but by choosing parameters from the posterior distributions, all these posterior quantities of interest can be approximated with arbitrary precision using the Monte Carlo method.
Let $g(θ)$ be any function. The law of large numbers says that if $θ_{(1)}, . . . , θ_{(S)}$ are $i.i.d.$ samples from $p(θ|y_1, . . . , y_n)$, then
$$ \frac{1}{S}\sum^S_{s=1}g(θ^{(s)}) → E[g(θ)|y_1, . . . , y_n] = \int g(θ)p(θ|y_1, . . . , y_n) dθ as S → ∞ . $$
This implies that as $S → ∞ ,$
The sequence {$(θ^{(1)}_1, θ^{(1)}_2), . . . ,(θ^{(S)}_1, θ^{(S)}_2)$} consists of S independent samples from the joint posterior distribution of $θ_1$ and $θ_2$, and can be used to make Monte Carlo approximations to posterior quantities of interest.