Posterior summaries from draws

Proposition taught in Bayesian foundations, part of Bayesian statistics

Posterior summaries from draws

Type Proposition
Taught in Bayesian foundations
Part Bayesian statistics
Worked examples 2
Questions using it 4

Let \theta^{(1)},\ldots,\theta^{(L)} be posterior draws and let \theta_{(1)}\leq\cdots\leq\theta_{(L)} be their order statistics. Then \begin{align*} \widehat{\mathbb{E}}(\theta\mid x) &=\frac1L\sum_{l=1}^L\theta^{(l)},\\ \widehat{\operatorname{Var}}(\theta\mid x) &=\frac1{L-1}\sum_{l=1}^L \left(\theta^{(l)}-\overline\theta\right)^2,\\ \widehat{\mathbb{P}}(\theta<c\mid x) &=\frac1L\sum_{l=1}^L\mathbf{1}\{\theta^{(l)}<c\}. \end{align*} The posterior median and equal-tailed interval are estimated by the corresponding empirical quantiles. With L=1000, the median can be estimated by the average of \theta_{(500)} and \theta_{(501)}, and a direct 95\% interval approximation by [\theta_{(25)},\theta_{(975)}]. The precise indices vary slightly with a software package’s quantile convention.

Examples

Q25. Gibbs, Metropolis–Hastings, and Monte Carlo error

Question page: Q25. Gibbs, Metropolis–Hastings, and Monte Carlo error

[25 marks]

Consider two random variables X and Y with joint density f(x,y)\propto x^3y\exp(-xy-2x-3y), \qquad x>0,\quad y>0, and zero density elsewhere.

Support clarification: the supplied screenshot does not state the support. This question uses the positive quadrant.

(a) [10 marks] Describe a Gibbs sampler to generate samples from the above joint density.

(b) [10 marks] Describe a random-walk Metropolis–Hastings algorithm to generate samples from the above joint density, using proposal increments distributed as N_2(\mathbf 0,\sigma^2I_2), where I_2 is the identity matrix and \sigma>0 is fixed.

(c) [5 marks] Briefly describe how you would obtain an approximate 95% confidence interval for \mathbb{E}(Y) from the simulation output.

Definitions. Take the joint support to be x,y>0. We use \operatorname{Gamma}(a,b) in shape–rate form, with kernel z^{a-1}e^{-bz} on z>0. N_2(\mathbf m,\Sigma) uses a covariance matrix \Sigma, and I_2 is the 2\times2 identity. The estimand in (c) is the fixed target mean \mu_Y=\mathbb{E}_f(Y); its simulation uncertainty is measured by a Monte Carlo standard error (MCSE).

(a) Gibbs sampler. For the conditional distribution of X, fix y>0 and retain all factors involving x: \begin{aligned} f(x\mid y)&\propto x^3y\exp(-xy-2x-3y)\\ &\propto x^3\exp\{-(y+2)x\},\qquad x>0. \end{aligned} The power 3=4-1 identifies shape 4, and the coefficient of x in the negative exponential is the rate y+2. Thus \boxed{X\mid Y=y\sim\operatorname{Gamma}(4,y+2).} Similarly, fixing x>0 gives f(y\mid x)\propto y\exp\{-(x+3)y\},\qquad y>0, and hence \boxed{Y\mid X=x\sim\operatorname{Gamma}(2,x+3).} The corresponding normalised conditional densities are f(x\mid y)=\frac{(y+2)^4}{\Gamma(4)}x^3e^{-(y+2)x}, \qquad f(y\mid x)=\frac{(x+3)^2}{\Gamma(2)}y e^{-(x+3)y}.

Choose any valid initial pair (x^{(0)},y^{(0)}) with both coordinates positive. A systematic-scan Gibbs sampler is:

1. For t=1,\ldots,T, draw X^{(t)}\sim\operatorname{Gamma}(4,Y^{(t-1)}+2).

2. Using this newly drawn value, draw Y^{(t)}\sim\operatorname{Gamma}(2,X^{(t)}+3).

3. Record the pair (X^{(t)},Y^{(t)}) and repeat.

Every draw is from an exact full conditional, so there is no accept–reject step. Each conditional update preserves the target joint distribution; their composition therefore has f(x,y) as an equilibrium distribution. After checking convergence and allowing suitable warm-up, retained pairs can be used to approximate target expectations. Successive pairs are generally dependent, not iid.

(b) Random-walk Metropolis–Hastings. Let the current state be (x,y) with x,y>0. Draw an increment \boldsymbol\varepsilon\sim N_2(\mathbf 0,\sigma^2I_2) and propose (x',y')=(x,y)+\boldsymbol\varepsilon. Equivalently, the proposed state is N_2((x,y)^\top,\sigma^2I_2). The zero-mean normal distribution describes the increment, not an independence proposal centred at the origin.

The proposal is symmetric: q\bigl((x',y')\mid(x,y)\bigr) =q\bigl((x,y)\mid(x',y')\bigr), because reversing the increment leaves its squared length unchanged. Thus the proposal densities cancel from the Metropolis–Hastings ratio. For a proposal in the positive quadrant, \begin{aligned} r&=\frac{f(x',y')}{f(x,y)}\\ &=\frac{(x')^3y' e^{-x'y'-2x'-3y'}}{x^3y e^{-xy-2x-3y}}\\ &=\left(\frac{x'}x\right)^3\frac{y'}y \exp\{-x'y'+xy-2(x'-x)-3(y'-y)\}. \end{aligned} Therefore the acceptance probability is \boxed{ a\bigl((x,y),(x',y')\bigr)= \begin{cases} \min\{1,r\},&x'>0,\ y'>0,\\ 0,&\text{otherwise}. \end{cases}} The unknown normalising constant of f cancels in the target-density ratio.

Start from a valid positive pair. For each iteration t=1,\ldots,T:

1. Draw an independent N_2(\mathbf0,\sigma^2I_2) increment and add it to (X^{(t-1)},Y^{(t-1)}) to obtain (x',y').

2. If either proposed coordinate is nonpositive, reject the proposal. Otherwise compute a above and draw U\sim U(0,1) independently; accept if U\leq a.

3. On acceptance set (X^{(t)},Y^{(t)})=(x',y'); on rejection set (X^{(t)},Y^{(t)})=(X^{(t-1)},Y^{(t-1)}).

4. Record the state, including repeated states caused by rejection.

An out-of-support proposal is rejected, not repeatedly redrawn until it is positive: such redrawing would change the proposal distribution and its acceptance ratio. Assess convergence before using the retained output.

(c) A 95% confidence interval for \mathbb{E}(Y). The target is the fixed number \mu_Y=\mathbb{E}_f(Y), estimated using a Monte Carlo average. After suitable warm-up and convergence checks, retain N consecutive values Y_1,\ldots,Y_N and calculate \bar Y_N=\frac1N\sum_{t=1}^{N}Y_t. These values are generally autocorrelated, so the iid standard error s_Y/\sqrt N is not generally valid.

One practical method is batch means. Write N=ab, with a non-overlapping batches of length b, and calculate \bar Y_k^{\mathrm{batch}}= \frac1b\sum_{t=(k-1)b+1}^{kb}Y_t,\qquad k=1,\ldots,a, s_B^2=\frac1{a-1}\sum_{k=1}^{a} \bigl(\bar Y_k^{\mathrm{batch}}-\bar Y_N\bigr)^2. With batches long relative to the chain’s dependence and sufficiently many batches, the batch averages are approximately independent and approximately normal under the usual Markov-chain central limit conditions. An estimated Monte Carlo standard error is \widehat{\operatorname{MCSE}}(\bar Y_N)=\frac{s_B}{\sqrt a}. An approximate 95% confidence interval is then \boxed{\bar Y_N\ \pm\ t_{0.975,a-1}\frac{s_B}{\sqrt a}.} For large a, the t quantile can be replaced by 1.96. The simulation should be long enough for the interval to be reasonably stable when the batch length is increased.

Equivalently, for a stationary chain satisfying a central limit theorem, \sqrt N(\bar Y_N-\mu_Y)\ \xrightarrow{d}\ N(0,\tau_Y^2), where the long-run variance is \tau_Y^2=\gamma_0+2\sum_{k=1}^{\infty}\gamma_k, \qquad \gamma_k=\operatorname{Cov}(Y_t,Y_{t+k}). Batch means estimates this variance by \widehat\tau_Y^2=b s_B^2, giving the same standard error \sqrt{\widehat\tau_Y^2/N}=s_B/\sqrt a.

Key result. This interval measures Monte Carlo error in estimating \mathbb{E}(Y). The empirical 2.5% and 97.5% quantiles of the Y_t values instead give an interval for the variable Y; they are not a confidence interval for its mean.

C11. Monte Carlo estimation

Question page: C11. Monte Carlo estimation

What is a Monte Carlo estimate of \mathbb{E}_\pi[h(\Theta)], and what is Monte Carlo error?

Definitions/notation. For a target distribution \pi and integrable function h, a Monte Carlo estimate is N^{-1}\sum_{t=1}^N h(\theta^{(t)}), using independent draws from \pi or suitable ergodic MCMC output. Monte Carlo error is the finite-simulation discrepancy from \mathbb{E}_\pi[h(\Theta)], not the posterior uncertainty about \Theta; its assessment must account for dependence when draws are correlated.

Used in

See also

Credible-region constructions · Bayesian credible interval · Bayes’ rule · Proper prior and proper posterior · Posterior point summaries · R’s distribution interface · Prior and posterior predictive distributions · Predictive-variance decomposition

Glossary and index