We may be interested in means and standard deviations of some function of $\theta$,. When comparing two or more populations we may be interested in the posterior distributionof $|\theta_{1} −\theta_{2}|, \theta_{1}/\theta_{1}$, or max$({\theta_{1}, . . . , \theta_{m}})$, all of which are functions of more
than one parameter. Obtaining exact values for these posterior quantities can
be difficult or impossible, but if we can generate random sample values of
the parameters from their posterior distributions, then all of these posterior
quantities of interest can be approximated to an arbitrary degree of precision
using the Monte Carlo method.
Let $\theta$ be a parameter of interest and let $y_{1}, . . . ,y_{n}$ be the numerical values of a sample from a distribution $p(y_{1}, . . . ,y_{n})$. Suppose we could sample some number $S$ of independent, random $\theta$-values from the posterior distribution $p(\theta|y_{1}, . . . ,y_{n})$:
$$\theta^{(1)}, . . . , \theta^{(s)} \sim p(\theta|y_{1}, . . . ,y_{n})$$.
Then the empirical distribution of the samples $\theta^{(1)}, . . . , \theta^{(s)}$ would approximate
$p(\theta|y_{1}, . . . ,y_{n})$, with the approximation improving with increasing $S$.The empirical distribution of $\theta^{(1)}, . . . , \theta^{(s)} $is known as a Monte Carlo approximation to $p(\theta|y_{1}, . . . ,y_{n})$.
Additionally, let $g(\theta)$ be (just about) any function of $\theta$. The law of large numbers says that if $\theta^{(1)}, . . . , \theta^{(s)} $ are i.i.d. samples from $p(\theta|y_{1}, . . . ,y_{n})$, then as $S\rightarrow \infty$
$$\frac{1}{S}\sum_{s=11}^{S}g(\theta) ^{(s)}\rightarrow E(g(\theta)|y_{1}, . . . ,y_{n}) $$
Example :
Suppose we are interested in the posterior distribution of some computable
function $g(\theta)$ of $\theta$. In the binomial model, for example, we are sometimes interested in the log odds :
$$log odds (\theta) = log \frac{\theta}{1-\theta}=\gamma$$
this quantity of interest can be computed using a Monte Carlo approach:
sample $\theta^{(1)} \sim p(\theta|y_{1}, . . . ,y_{n})$, compute
$\gamma^{(1)} = g(\theta^{(1)})$
sample $\theta^{(2)} \sim p(\theta|y_{1}, . . . ,y_{n})$, compute
$\gamma^{(2)} = g(\theta^{(2)})$
.
.
.
sample $\theta^{(s)} \sim p(\theta|y_{1}, . . . ,y_{n})$, compute
$\gamma^{(s)} = g(\theta^{(s)})$
Then the sequence $\{\gamma^{(s)}, . . . , \gamma^{(s)}\}$ constitutes $S$ independent samples from $p(\gamma|y_{1}, . . . ,y_{n})$.
Let $\theta$ be a population proportion . Using a binomial sampling model and a uniform prior distribution the posterior distribution of $\theta$ will be beta distributed say with beta(45, 420). Using the Monte Carlo algorithm described above, we can obtain samples of the log-odds
from the posterior distribution of $\theta$ . In R, the Monte Carlo algorithm involves only a few commands:
set.seed(1)
a<-1 ; b<-1 # Parameter for the prior
# Prior for theta
theta.prior.sim<-rbeta(10000,a,b) # Beta(1,1) is also uniform
# Prior of the function of interest log odds (gamma)
gamma.prior.sim<- log( theta.prior.sim/(1-theta.prior.sim) )
# Posterior for theta
theta.post.sim<-rbeta(10000,45,420)
# Posterior for log odds (gamma)
gamma.post.sim<- log( theta.post.sim/(1-theta.post.sim) )
# Plot the prior and the posterior for gammma
par(mar=c(3,3,1,1),mgp=c(1.75,.75,0))
par(mfrow=c(2,3))
par(cex=.8)
par(mfrow=c(1,2),mar=c(3,3,1,1), mgp=c(1.75,.75,.0))
plot(density(gamma.prior.sim,adj=2),xlim=c(-5,5),main="",
xlab=expression(gamma),
ylab=expression(italic(p(gamma))),col="gray")
plot(density(gamma.post.sim,adj=2),xlim=c(-5,5),main="",xlab=expression(gamma),
ylab=expression(paste(italic("p("),gamma,"|",y[1],"...",y[n],")",
sep="")) )
lines(density(gamma.prior.sim,adj=2),col="gray")
For further reading I recommend A First Course in Bayesian Statistical Methods
and Doing Bayesian Data Analysis: A Tutorial Introduction with R and BUGS
