Bayesian Inference
This page condenses UNSW MATH5960 (Bayesian Inference and Computation). It assumes the Probability page and sits under Statistics, whose table of distributions it uses throughout. Gamma distributions are in the rate parameterisation: \(\mathrm{Gamma}(a,b)\) has density \(\frac{b^a}{\Gamma(a)}\theta^{a-1}e^{-b\theta}\), mean \(a/b\) and variance \(a/b^2\).
The course has two halves that keep meeting. One is modelling: write a likelihood and a prior, get a posterior, and answer every question with an integral against it. The other is computation: when those integrals have no closed form, replace them by averages over draws, first independent (Monte Carlo, inversion, rejection, importance sampling), then dependent (Gibbs and Metropolis-Hastings). The last third uses both halves on model comparison, hierarchical models and mixtures.
The Bayesian Recipe
Frequentist inference treats \(\theta\) as fixed and the data as random; its statements (a confidence interval covers \(\theta\) in 95% of repeated experiments) are about procedures. Bayesian inference treats \(\theta\) as uncertain, describes that uncertainty with a distribution, and conditions on the data actually observed.
For events \(A\) and a partition \(B_1,\dots,B_k\) of the sample space with \(P(B_j)>0\), \[ P(B_j\mid A)=\frac{P(A\mid B_j)P(B_j)}{\sum_{i=1}^kP(A\mid B_i)P(B_i)} . \] For densities, \(\pi(\theta\mid x)=f(x\mid\theta)\pi(\theta)\big/\int f(x\mid\vartheta)\pi(\vartheta)\,d\vartheta\).
\(P(B_j\mid A)P(A)=P(A\cap B_j)=P(A\mid B_j)P(B_j)\), and the law of total probability gives \(P(A)=\sum_iP(A\mid B_i)P(B_i)\) because the \(A\cap B_i\) are disjoint with union \(A\). The density version is the same computation with the sum replaced by an integral.
Given a model \(L(x\mid\theta)\) for data \(x\) (the likelihood, read as a function of \(\theta\)) and a prior density \(\pi(\theta)\) on the parameter space \(\Theta\), the marginal likelihood (normalising constant, evidence) is \[ m(x)=\int_\Theta L(x\mid\theta)\pi(\theta)\,d\theta , \] and the posterior is \[ \pi(\theta\mid x)=\frac{L(x\mid\theta)\pi(\theta)}{m(x)}\ \propto\ L(x\mid\theta)\pi(\theta). \]
Almost every calculation on this page uses the proportional form. Write prior times likelihood, drop every factor that does not involve \(\theta\), and recognise the kernel of a known density; the constant then comes for free. Only the marginal likelihood itself (Bayes factors, below) needs all the constants kept.
- Specify the likelihood \(L(x\mid\theta)\).
- Specify the prior \(\pi(\theta)\).
- Compute the posterior \(\pi(\theta\mid x)\propto L(x\mid\theta)\pi(\theta)\).
- Draw inferences from the posterior: point estimates, intervals, probabilities of hypotheses, predictions.
A set \(C\subseteq\Theta\) is a \(100(1-\alpha)\%\) credible region if \(P(\theta\in C\mid x)=1-\alpha\). The equal-tailed interval runs between the \(\alpha/2\) and \(1-\alpha/2\) posterior quantiles. The high density region (HDR) is \(C=\{\theta:\pi(\theta\mid x)\ge k\}\) with \(k\) chosen so the posterior mass is \(1-\alpha\); for a unimodal posterior it is the shortest credible interval.
Unlike a confidence interval, a credible interval means what people usually want it to mean: given these data, \(\theta\) lies in \(C\) with probability \(1-\alpha\). For a skewed posterior the equal-tailed interval can exclude the most probable values, so the HDR is the one to report.
For a future observation \(y\), conditionally independent of \(x\) given \(\theta\), \[ p(y\mid x)=\int_\Theta f(y\mid\theta)\,\pi(\theta\mid x)\,d\theta . \] The prior predictive \(p(y)=\int f(y\mid\theta)\pi(\theta)\,d\theta\) is the same average taken under the prior; comparing it with the data checks whether the prior is sensible on the scale of observables.
A hospital has \(r=0\) deaths in \(n=47\) operations, \(r\sim\mathrm{Bin}(n,\theta)\). The likelihood \((1-\theta)^{47}\) is strictly decreasing, so \(\hat\theta=0\) sits on the boundary; the plug-in variance \(\hat\theta(1-\hat\theta)/n\) is \(0\), the Wald interval is \([0,0]\), and the plug-in predictive says the next baby survives with probability one.
With a \(\mathrm{Unif}(0,1)\) prior the posterior is \(\pi(\theta\mid r)\propto(1-\theta)^{47}\), that is \(\mathrm{Beta}(1,48)\). It is decreasing, so the 95% HDR is \([0,c]\) with \(1-(1-c)^{48}=0.95\), giving \(c=1-0.05^{1/48}\approx0.0605\). The posterior probability that \(\theta\le0.05\) is \(1-0.95^{48}\approx0.915\), and the predictive probability that the next operation is a death is \(E(\theta\mid r)=1/49\approx0.020\): small but not zero. No asymptotics were needed.
Conjugate Priors
A family \(\mathcal P\) of priors is conjugate for a likelihood if every prior in \(\mathcal P\) gives a posterior in \(\mathcal P\). Updating is then arithmetic on the family’s parameters.
If \(x\mid\theta\sim\mathrm{Bin}(n,\theta)\) and \(\theta\sim\mathrm{Beta}(p,q)\), then \(\theta\mid x\sim\mathrm{Beta}(p+x,\,q+n-x)\).
\(\pi(\theta\mid x)\propto\theta^x(1-\theta)^{n-x}\cdot\theta^{p-1}(1-\theta)^{q-1}=\theta^{p+x-1}(1-\theta)^{q+n-x-1}\), the Beta kernel.
If \(x_1,\dots,x_n\mid\theta\) are iid \(\mathrm{Pois}(\theta)\) and \(\theta\sim\mathrm{Gamma}(\alpha,\beta)\), then \[ \theta\mid x\sim\mathrm{Gamma}(a,b),\qquad a=\alpha+\textstyle\sum_ix_i,\quad b=\beta+n, \] and the predictive distribution of a new count is negative binomial, \[ p(y\mid x)=\binom{a+y-1}{y}\Big(\frac{b}{b+1}\Big)^a\Big(\frac{1}{b+1}\Big)^y,\qquad y=0,1,2,\dots \]
The posterior: \(\prod_i\theta^{x_i}e^{-\theta}\cdot\theta^{\alpha-1}e^{-\beta\theta}=\theta^{\alpha+\sum x_i-1}e^{-(\beta+n)\theta}\). The predictive uses the Gamma integral \(\int_0^\infty\theta^{s-1}e^{-b\theta}\,d\theta=\Gamma(s)/b^s\): \[ p(y\mid x)=\int_0^\infty\frac{\theta^ye^{-\theta}}{y!}\,\frac{b^a\theta^{a-1}e^{-b\theta}}{\Gamma(a)}\,d\theta =\frac{b^a}{y!\,\Gamma(a)}\cdot\frac{\Gamma(a+y)}{(b+1)^{a+y}}, \] and \(\Gamma(a+y)/(y!\,\Gamma(a))=\binom{a+y-1}{y}\).
The posterior mean \(a/b=\frac{\alpha+\sum x_i}{\beta+n}\) is a weighted average of the prior mean \(\alpha/\beta\) (weight \(\beta/(\beta+n)\)) and \(\bar x\) (weight \(n/(\beta+n)\)): the prior acts like \(\beta\) extra observations with total \(\alpha\). The same reading works for every conjugate pair. In the notes the pair is applied to 65 years of Melbourne per-capita water use; the analysis is easy and plainly wrong, because a constant Poisson rate ignores the trend, which is a reminder that conjugacy is a computational convenience and not a modelling argument.
Let \(x_1,\dots,x_n\mid\theta\) be iid \(N(\theta,\sigma^2)\) with data precision \(\tau=1/\sigma^2\) known, and \(\theta\sim N(b,1/c)\). Then \[ \theta\mid x\sim N\Big(\frac{cb+n\tau\bar x}{c+n\tau},\ \frac{1}{c+n\tau}\Big). \] Equivalently, posterior precision is prior precision plus \(n\) times data precision, and the posterior mean is \(\gamma_nb+(1-\gamma_n)\bar x\) with \(\gamma_n=c/(c+n\tau)\). The posterior depends on the data only through \(\bar x\) (sufficiency), and both \(n\to\infty\) and \(c\to0\) give \(\theta\mid x\approx N(\bar x,\sigma^2/n)\).
Up to factors free of \(\theta\), \[ \pi(\theta\mid x)\propto\exp\Big(-\frac c2(\theta-b)^2-\frac\tau2\sum_i(x_i-\theta)^2\Big) \propto\exp\Big(-\frac12\big[A\theta^2-2B\theta\big]\Big), \] with \(A=c+n\tau\) and \(B=cb+n\tau\bar x\). Completing the square, \(A\theta^2-2B\theta=A(\theta-B/A)^2-B^2/A\), and the last term is free of \(\theta\). So the posterior is normal with mean \(B/A\) and variance \(1/A\).
Suppose \(f(x\mid\theta)=h(x)\,g(\theta)\exp\{c(\theta)t(x)\}\). Then the family \[ \pi(\theta)\propto g(\theta)^d\exp\{b\,c(\theta)\} \] is conjugate, and \(n\) observations update \(d\mapsto d+n\) and \(b\mapsto b+\sum_it(x_i)\).
\(L(x\mid\theta)\pi(\theta)\propto\prod_ih(x_i)\,g(\theta)^n e^{c(\theta)\sum_it(x_i)}\cdot g(\theta)^de^{bc(\theta)}\propto g(\theta)^{d+n}e^{(b+\sum t(x_i))c(\theta)}\).
The standard pairs, all instances of the proposition or a two-line calculation like the Beta-binomial one:
| likelihood | prior | posterior |
|---|---|---|
| \(x\sim\mathrm{Bin}(n,\theta)\) | \(\mathrm{Beta}(p,q)\) | \(\mathrm{Beta}(p+x,\,q+n-x)\) |
| \(x_i\sim\mathrm{Pois}(\theta)\) | \(\mathrm{Gamma}(p,q)\) | \(\mathrm{Gamma}(p+\sum x_i,\,q+n)\) |
| \(x_i\sim N(\theta,1/\tau)\), \(\tau\) known | \(N(b,1/c)\) | \(N\big(\frac{cb+n\tau\bar x}{c+n\tau},\frac1{c+n\tau}\big)\) |
| \(x_i\sim\mathrm{Gamma}(k,\theta)\), \(k\) known | \(\mathrm{Gamma}(p,q)\) | \(\mathrm{Gamma}(p+nk,\,q+\sum x_i)\) |
| \(x_i\sim\mathrm{Geom}(\theta)\) (success prob.) | \(\mathrm{Beta}(p,q)\) | \(\mathrm{Beta}(p+n,\,q+\sum x_i-n)\) |
| \(x\sim\mathrm{NegBin}(r,\theta)\) | \(\mathrm{Beta}(p,q)\) | \(\mathrm{Beta}(p+r,\,q+x)\) |
Let \(\pi_1,\dots,\pi_k\) be conjugate priors with posteriors \(\pi_i(\theta\mid x)\) and marginal likelihoods \(m_i(x)=\int L(x\mid\theta)\pi_i(\theta)\,d\theta\). Under the mixture prior \(\pi=\sum_iw_i\pi_i\), \[ \pi(\theta\mid x)=\sum_{i=1}^kw_i^*\,\pi_i(\theta\mid x),\qquad w_i^*=\frac{w_i\,m_i(x)}{\sum_jw_j\,m_j(x)} . \]
\(\pi(\theta\mid x)=\sum_iw_iL(x\mid\theta)\pi_i(\theta)/m(x)\) with \(m(x)=\sum_jw_jm_j(x)\). Write \(L\pi_i=m_i(x)\pi_i(\theta\mid x)\) in each term.
The components update as they would alone, and each weight is rescaled by how well its component predicted the data. Mixtures make the conjugate family flexible enough to express a bimodal belief or a sceptical prior with a heavy tail, at no computational cost.
Choosing a Prior
A prior with \(\int_\Theta\pi(\theta)\,d\theta=\infty\) is improper. It may be used only when the posterior it produces is proper, \(\int L(x\mid\theta)\pi(\theta)\,d\theta<\infty\), which can depend on the observed data.
The maximum a posteriori estimate \(\hat\theta_{\mathrm{MAP}}=\arg\max_\theta\{\log L(x\mid\theta)+\log\pi(\theta)\}\) equals the MLE when \(\pi\) is flat on a set containing the MLE.
A flat prior is not a statement of ignorance, because flatness is not preserved by a change of variable: if \(\theta\sim\mathrm{Unif}(0,1)\) then \(\phi=\theta^2\) has density \(1/(2\sqrt\phi)\) and \(P(\phi\le\tfrac14)=\tfrac12\). Jeffreys asked for a rule that gives the same answer whichever parameterisation one starts from.
\[ \pi_J(\theta)\propto|I(\theta)|^{1/2},\qquad I(\theta)=-E\Big[\frac{\partial^2}{\partial\theta^2}\log L(x\mid\theta)\Big], \] the expectation taken over \(x\sim L(\cdot\mid\theta)\); for vector \(\theta\), \(I\) is a matrix and \(|\cdot|\) its determinant.
If \(\phi=g(\theta)\) is monotone and differentiable, then \(I(\phi)=I(\theta)(d\theta/d\phi)^2\), so transforming \(\pi_J(\theta)\) as a density gives Jeffreys’ prior for \(\phi\).
Let \(\ell(\theta)=\log L(x\mid\theta)\). By the chain and product rules, \[ \frac{d^2\ell}{d\phi^2}=\ell’’(\theta)\Big(\frac{d\theta}{d\phi}\Big)^2+\ell’(\theta)\frac{d^2\theta}{d\phi^2}. \] Take expectations over \(x\). The score has mean zero, \(E[\ell’(\theta)]=\int\partial_\theta L\,dx=\partial_\theta\int L\,dx=0\), so the second term vanishes and \(I(\phi)=I(\theta)(d\theta/d\phi)^2\). Hence \(\pi_J(\theta)\,|d\theta/d\phi|\propto I(\theta)^{1/2}|d\theta/d\phi|=I(\phi)^{1/2}\).
- Binomial: \(\ell=x\log\theta+(n-x)\log(1-\theta)\), \(I(\theta)=\frac{n}{\theta(1-\theta)}\), so \(\pi_J=\mathrm{Beta}(\tfrac12,\tfrac12)\), which is proper.
- Normal mean, \(\sigma^2\) known: \(I\) is constant, so \(\pi_J\propto1\).
- Normal variance \(\theta=\sigma^2\), mean known: \(I(\theta)=\frac{n}{2\theta^2}\), so \(\pi_J\propto1/\theta\).
- Normal mean and variance jointly: the determinant of the \(2\times2\) information gives \(\pi_J\propto\sigma^{-3}\), which disagrees with the product \(1\times1/\sigma^2\) of the two one-parameter answers.
The objections: \(\pi_J\) depends on the sampling model and not just on the parameter (so it violates the likelihood principle), it is often improper, and the joint and separate versions disagree in the normal model. In practice one uses Jeffreys’ prior or a weakly informative proper prior and checks that the posterior is proper and not sensitive to the choice. The prior predictive check is the practical test: simulate data from \(p(y)\) and ask whether they look like anything the world could produce.
Multiparameter Models
For \(\theta=(\theta_1,\theta_2)\), the joint posterior is \(\pi(\theta_1,\theta_2\mid x)\propto L(x\mid\theta_1,\theta_2)\pi(\theta_1,\theta_2)\) and the marginal posterior of the parameter of interest is \[ \pi(\theta_1\mid x)=\int\pi(\theta_1,\theta_2\mid x)\,d\theta_2=\int\pi(\theta_1\mid\theta_2,x)\,\pi(\theta_2\mid x)\,d\theta_2 . \] The nuisance parameter \(\theta_2\) is averaged over rather than plugged in. With posterior draws, marginalising means discarding the other coordinates.
Let \(y_1,\dots,y_n\) be iid \(N(\mu,\sigma^2)\) with the prior \(\pi(\mu,\sigma^2)\propto1/\sigma^2\), and \(s^2=\frac1{n-1}\sum(y_i-\bar y)^2\). Then \[ \mu\mid\sigma^2,y\sim N(\bar y,\sigma^2/n),\qquad \sigma^2\mid y\sim\mathrm{InvGamma}\Big(\frac{n-1}2,\frac{(n-1)s^2}2\Big),\qquad \frac{\mu-\bar y}{s/\sqrt n}\,\Big|\,y\sim t_{n-1}. \]
The decomposition \(\sum(y_i-\mu)^2=(n-1)s^2+n(\bar y-\mu)^2\) gives \[ \pi(\mu,\sigma^2\mid y)\propto(\sigma^2)^{-n/2-1}\exp\Big(-\frac{(n-1)s^2+n(\bar y-\mu)^2}{2\sigma^2}\Big). \] As a function of \(\mu\) this is the \(N(\bar y,\sigma^2/n)\) kernel. Integrating it out contributes \(\sqrt{2\pi\sigma^2/n}\), leaving \((\sigma^2)^{-(n-1)/2-1}e^{-(n-1)s^2/(2\sigma^2)}\), the inverse Gamma kernel. Integrating \(\sigma^2\) out instead with the Gamma integral gives \(\pi(\mu\mid y)\propto[(n-1)s^2+n(\mu-\bar y)^2]^{-n/2}\propto\big[1+\frac{1}{n-1}\big(\frac{\mu-\bar y}{s/\sqrt n}\big)^2\big]^{-n/2}\), the \(t_{n-1}\) kernel.
The Bayesian 95% interval for \(\mu\) is then \(\bar y\pm t_{n-1,0.975}\,s/\sqrt n\), numerically the classical \(t\) interval, with the Bayesian reading. The fully conjugate prior \(\mu\mid\sigma^2\sim N(\mu_0,\sigma^2/\kappa_0)\), \(\sigma^2\sim\mathrm{InvGamma}(\alpha_0,\beta_0)\) updates to \[ \kappa_n=\kappa_0+n,\quad \mu_n=\frac{\kappa_0\mu_0+n\bar y}{\kappa_n},\quad \alpha_n=\alpha_0+\frac n2,\quad \beta_n=\beta_0+\frac{(n-1)s^2}{2}+\frac{\kappa_0n(\bar y-\mu_0)^2}{2\kappa_n}. \]
For \(y_1,\dots,y_n\) iid \(N_p(\mu,\Sigma)\):
- with \(\Sigma\) known and \(\mu\sim N_p(\mu_0,\Lambda_0)\), the posterior is \(N_p(\mu_n,\Lambda_n)\) with \(\Lambda_n^{-1}=\Lambda_0^{-1}+n\Sigma^{-1}\) and \(\mu_n=\Lambda_n(\Lambda_0^{-1}\mu_0+n\Sigma^{-1}\bar y)\), the matrix form of “precisions add”;
- with \(\mu\) known, an inverse Wishart prior on \(\Sigma\) is conjugate, \(\mathrm{IW}(\nu_0,S_0)\mapsto\mathrm{IW}\big(\nu_0+n,\ S_0+\sum_i(y_i-\mu)(y_i-\mu)^\top\big)\), and correspondingly a Wishart prior on the precision \(\Sigma^{-1}\).
The Multivariate Analysis page treats the Wishart distribution itself.
Decisions and Point Estimates
A loss function \(L(\theta,d)\) is the cost of taking decision \(d\) when the parameter is \(\theta\). The Bayes decision minimises the expected posterior loss \[ \rho(d\mid x)=E[L(\theta,d)\mid x]=\int_\Theta L(\theta,d)\,\pi(\theta\mid x)\,d\theta . \]
Point estimation is the case \(d\in\Theta\), and the four standard losses produce the four standard summaries.
Under \(L(\theta,d)=(\theta-d)^2\) the Bayes estimate is \(E(\theta\mid x)\).
With \(\bar\theta=E(\theta\mid x)\), \(E[(\theta-d)^2\mid x]=\mathrm{Var}(\theta\mid x)+(\bar\theta-d)^2\), since the cross term \(2(\bar\theta-d)E[\theta-\bar\theta\mid x]\) is zero. Minimise the second term.
Under \(L(\theta,d)=g(d-\theta)\) for \(d>\theta\) and \(h(\theta-d)\) for \(d<\theta\) (cost \(g\) per unit of overestimate, \(h\) per unit of underestimate), the Bayes estimate is the \(\frac{h}{g+h}\) quantile of the posterior. With \(g=h\) (absolute error loss) it is the posterior median.
Let \(F\) be the posterior cdf. Then \[ \rho(d)=g\int_{-\infty}^d(d-\theta)\,dF(\theta)+h\int_d^{\infty}(\theta-d)\,dF(\theta), \] and differentiating under the integrals (the boundary terms vanish because the integrands are zero at \(\theta=d\)) gives \(\rho’(d)=gF(d)-h(1-F(d))\). This is increasing in \(d\), so \(\rho\) is convex and is minimised where \(F(d)=h/(g+h)\).
Under \(L(\theta,d)=\mathbf 1\{|\theta-d|>\epsilon\}\), \(\rho(d)=1-P(|\theta-d|\le\epsilon\mid x)\approx1-2\epsilon\,\pi(d\mid x)\) for small \(\epsilon\); as \(\epsilon\to0\) the Bayes estimate is the posterior mode.
Decisions need not be estimates. A loaf costs \(c\) to bake and sells for \(s>c\); demand \(D\) is uncertain with posterior \(\pi(D\mid x)\). Measured against the ideal of baking exactly \(D\), baking \(b\) loaves loses \(c(b-D)\) when \(b>D\) (wasted loaves) and \((s-c)(D-b)\) when \(b\le D\) (missed margin). This is linear loss with \(g=c\) and \(h=s-c\), so the optimal batch is the \(\frac{s-c}{s}\) posterior quantile of demand. With \(c=1\), \(s=4\) and demand uniform on \(\{0,\dots,9\}\), the target level is \(\frac34\) and \(b^*=7\), well above the mean demand of \(4.5\): a missed sale costs three times a wasted loaf.
Posterior Asymptotics
Suppose the model is correct, with true value \(\theta_0\), identifiable, and \(\pi\) positive in a neighbourhood of \(\theta_0\). As \(n\to\infty\):
- (consistency) the posterior mass of every neighbourhood of \(\theta_0\) tends to \(1\);
- (Bernstein-von Mises) \(\pi(\theta\mid x)\) is approximately \(N\big(\hat\theta,\,I_n(\hat\theta)^{-1}\big)\), where \(\hat\theta\) is the MLE (or posterior mode) and \(I_n\) the observed information.
If \(x\sim f(\cdot\mid\theta_0)\) then \(\theta_0\) maximises \(\theta\mapsto E_{\theta_0}[\log f(x\mid\theta)]\).
By Jensen’s inequality for the concave \(\log\), \[ E_{\theta_0}\Big[\log\frac{f(x\mid\theta)}{f(x\mid\theta_0)}\Big]\le\log E_{\theta_0}\Big[\frac{f(x\mid\theta)}{f(x\mid\theta_0)}\Big]=\log\int f(x\mid\theta)\,dx=\log1=0 . \]
Why the theorem holds (sketch). Write \(\pi(\theta\mid x)\propto\pi(\theta)\exp\{n\ell_n(\theta)\}\) with \(\ell_n(\theta)=\frac1n\sum_i\log f(x_i\mid\theta)\). By the law of large numbers \(\ell_n(\theta)\to E_{\theta_0}[\log f(x\mid\theta)]\), which the lemma says peaks uniquely at \(\theta_0\) when the model is identifiable; the factor \(e^{n\ell_n}\) then crushes all mass away from \(\theta_0\), and the fixed prior cannot resist unless \(\pi(\theta_0)=0\). For normality, expand the log posterior to second order about its mode: the linear term vanishes there, and the quadratic coefficient is \(-\frac12I_n(\hat\theta)\) up to the \(O(1)\) prior contribution, which is negligible against the \(O(n)\) likelihood curvature. So every point estimate of the previous section converges to the MLE, and a credible interval and a confidence interval agree to first order.
Monte Carlo
Every posterior summary is an integral, \(E[g(\theta)\mid x]=\int g(\theta)\pi(\theta\mid x)\,d\theta\): means, variances, probabilities (\(g\) an indicator), predictive densities. Monte Carlo replaces the integral by an average over draws, and the rest of the course is about where the draws come from.
If \(\theta^{(1)},\dots,\theta^{(N)}\) are iid draws from \(\pi\), estimate \(\mu=E_\pi[g(\theta)]\) by \(\hat\mu_N=\frac1N\sum_{i=1}^Ng(\theta^{(i)})\). To integrate \(\int_a^bh(x)\,dx\) against no density at all, write it as \((b-a)E[h(U)]\) with \(U\sim\mathrm{Unif}(a,b)\); an infinite range is first mapped to \([0,1]\), for instance by \(x=u/(1-u)\).
If \(E_\pi|g(\theta)|<\infty\), the strong law of large numbers gives \(\hat\mu_N\to\mu\) almost surely. If also \(\sigma_g^2=\mathrm{Var}_\pi(g(\theta))<\infty\), then \(\mathrm{Var}(\hat\mu_N)=\sigma_g^2/N\) and, by the central limit theorem, \(\hat\mu_N\approx N(\mu,\sigma_g^2/N)\).
The Monte Carlo standard error \(s_g/\sqrt N\), with \(s_g\) the sample standard deviation of the \(g(\theta^{(i)})\), is reported alongside every estimate. For error at most \(E\) with probability \(1-\alpha\), take \(N\ge(z_{1-\alpha/2}\,s_g/E)^2\). The rate \(1/\sqrt N\) does not depend on the dimension of \(\theta\), which is why Monte Carlo beats quadrature in high dimensions, and the same draws serve every \(g\) at once: a marginal is obtained by discarding coordinates, and a predictive draw by drawing \(y^{(i)}\sim f(y\mid\theta^{(i)})\).
Inversion and rejection
If \(F\) is a continuous, strictly increasing cdf and \(U\sim\mathrm{Unif}(0,1)\), then \(X=F^{-1}(U)\) has cdf \(F\). (With the generalised inverse \(F^{-1}(u)=\inf\{x:F(x)\ge u\}\) it holds for every cdf.)
\(P(X\le x)=P(F^{-1}(U)\le x)=P(U\le F(x))=F(x)\), using that \(F\) is increasing.
For \(\mathrm{Exp}(\lambda)\), \(F(x)=1-e^{-\lambda x}\) and \(F^{-1}(u)=-\log(1-u)/\lambda\). Since \(1-U\) is also uniform, \(-\log(U)/\lambda\) will do.
Let \(f\) be a target density, known perhaps only up to a constant, and \(g\) a proposal density we can sample, with \(K\ge\sup_xf(x)/g(x)<\infty\). Repeat: draw \(x\sim g\) and \(u\sim\mathrm{Unif}(0,1)\); accept \(x\) if \(u\le f(x)/(Kg(x))\).
Accepted draws have density \(f/\int f\), and each proposal is accepted with probability \(\int f(x)\,dx/K\).
\[ P(X\le t,\ \text{accept})=\int_{-\infty}^tg(x)\,\frac{f(x)}{Kg(x)}\,dx=\frac1K\int_{-\infty}^tf(x)\,dx . \] Letting \(t\to\infty\) gives \(P(\text{accept})=\int f/K\), and dividing gives \(P(X\le t\mid\text{accept})=\int_{-\infty}^tf\big/\int f\).
The efficiency is the area under \(f\) over the area under the envelope \(Kg\), so the envelope should hug the target. In high dimensions no envelope does, and the acceptance rate collapses; that is the motivation for importance sampling and for MCMC.
Importance sampling
Let \(g\) be a density with \(g>0\) wherever \(f\,h\neq0\), draw \(x^{(i)}\sim g\) iid, and set weights \(w^{(i)}=f(x^{(i)})/g(x^{(i)})\).
- If \(f\) is normalised, \(\frac1N\sum_iw^{(i)}h(x^{(i)})\) is unbiased for \(E_f[h(X)]\).
- If only \(\tilde f\propto f\) is known, the self-normalised estimator \(\sum_iW^{(i)}h(x^{(i)})\) with \(W^{(i)}=\tilde w^{(i)}/\sum_j\tilde w^{(j)}\) is consistent but biased.
- \(\frac1N\sum_i\tilde f(x^{(i)})/g(x^{(i)})\) is unbiased for the normalising constant \(Z=\int\tilde f\).
(1) and (3): \(E_g\big[\frac{f(X)}{g(X)}h(X)\big]=\int\frac{f}{g}h\,g\,dx=\int fh\,dx\), with \(h\equiv1\) and \(\tilde f\) in place of \(f\) for (3). (2): the estimator is the ratio \(\frac1N\sum\tilde w^{(i)}h(x^{(i)})\big/\frac1N\sum\tilde w^{(i)}\); numerator and denominator converge almost surely to \(Z\,E_f[h]\) and \(Z\), so the ratio converges to \(E_f[h]\). A ratio of unbiased estimators is not unbiased in general.
With normalised weights \(W^{(i)}\), \(\mathrm{ESS}=1\big/\sum_i(W^{(i)})^2\). It lies in \([1,N]\): it equals \(N\) when all weights are equal and \(1\) when one weight carries everything. As \(N\to\infty\), \(\mathrm{ESS}/N\to1/E_g[w^2]\) for normalised \(w=f/g\).
The \(W^{(i)}\) are non-negative and sum to one. So \(\sum W_i^2\le(\sum W_i)^2=1\), and by Cauchy-Schwarz \(1=(\sum W_i\cdot1)^2\le N\sum W_i^2\). Invert.
A good proposal resembles the target in shape and has heavier tails. A proposal with lighter tails than \(f\) can give \(E_g[w^2]=\infty\): the estimator is still consistent, but its variance is infinite and the ESS is dominated by rare huge weights.
Markov Chain Monte Carlo
MCMC gives up independence: it builds a Markov chain whose stationary distribution is the posterior and uses one long path of the chain as a dependent sample. The theory needed is exactly what licenses path averages. For the probability of Markov chains in general see Markov chains and hidden Markov models.
A Markov chain on a state space \(E\) has \(P(X_{t+1}\in A\mid X_t=x,X_{t-1},\dots)=P(x,A)\), a transition kernel \(P\); for a finite chain the \(n\)-step probabilities are the entries of the matrix power \(P^n\). A distribution \(\pi\) is stationary (invariant) if \(X_t\sim\pi\) implies \(X_{t+1}\sim\pi\): \(\pi P=\pi\), or \(\int\pi(x)P(x,y)\,dx=\pi(y)\).
A chain is irreducible (\(\pi\)-irreducible on a general space) if from every state it can reach every state (every set of positive \(\pi\)-mass) with positive probability. The period of \(x\) is \(\gcd\{n\ge1:P^n(x,x)>0\}\); a chain is aperiodic if every period is \(1\). A self-loop \(P(x,x)>0\) makes \(x\) aperiodic.
If the chain is \(\pi\)-irreducible with stationary distribution \(\pi\) and \(\int|g|\,d\pi<\infty\), then \(\frac1N\sum_{t=1}^Ng(X_t)\to\int g\,d\pi\) for \(\pi\)-almost every starting value. If it is also aperiodic, then \(\|P^n(X_0,\cdot)-\pi\|_{TV}=\sup_A|P^n(X_0,A)-\pi(A)|\to0\): the chain’s distribution converges to \(\pi\) from almost any start.
The ergodic theorem needs irreducibility but not aperiodicity: a deterministic 3-cycle has time averages that converge though \(X_t\) never settles. So the recipe is: design a chain with the posterior as stationary distribution, check irreducibility (and aperiodicity if you want the draws themselves to be distributed as \(\pi\)), start anywhere, and average.
Weather is rain or sun with \(P(R,R)=0.4\), \(P(R,S)=0.6\), \(P(S,R)=0.2\), \(P(S,S)=0.8\). Solving \(\pi P=\pi\) with \(\pi( R)+\pi(S)=1\) gives \(0.6\,\pi( R)=0.2\,\pi(S)\), so \(\pi=(\tfrac14,\tfrac34)\). The second eigenvalue of \(P\) is \(0.4-0.2=0.2\), and from \(X_0=R\) the total variation distance is \(\tfrac34(0.2)^n\): geometric convergence, at rate given by the second eigenvalue.
The Gibbs sampler
For \(\theta=(\theta_1,\dots,\theta_p)\) with target \(\pi(\theta)\), one sweep replaces each component in turn by a draw from its full conditional given the newest values of the others: \[ \theta_j^{(t+1)}\sim\pi\big(\theta_j\mid\theta_1^{(t+1)},\dots,\theta_{j-1}^{(t+1)},\theta_{j+1}^{(t)},\dots,\theta_p^{(t)}\big),\qquad j=1,\dots,p . \] Variants update the components in random order, or one random component per step.
To find a full conditional, write the joint posterior, treat the other components as constants, drop every factor free of \(\theta_j\), and recognise the kernel. This is where conjugacy pays even in models that are not conjugate as a whole.
For \(y_i\) iid \(N(\mu,\sigma^2)\) with prior \(1/\sigma^2\), the joint posterior above gives the full conditionals \[ \mu\mid\sigma^2,y\sim N(\bar y,\sigma^2/n),\qquad \sigma^2\mid\mu,y\sim\mathrm{InvGamma}\Big(\frac n2,\ \frac12\sum_i(y_i-\mu)^2\Big), \] and alternating them samples the joint posterior, though here it could also be sampled exactly by composition.
If \(\pi(\theta)>0\) whenever every marginal \(\pi(\theta_j)\) is positive (the positivity condition), the full conditionals determine the joint distribution uniquely. The condition is sufficient, not necessary.
Let \(\pi\) be uniform on the union of the two squares \([0,1]^2\) and \([2,3]^2\). Each axis-parallel Gibbs move stays inside the current square, so the chain is reducible and its limit depends on the start; positivity fails, since both marginals are positive on \([2,3]\times[0,1]\) where the joint is zero. A milder failure is correlation: for a bivariate normal with correlation \(\rho\), the Gibbs chain for one coordinate is an autoregression with coefficient \(\rho^2\) per sweep, so mixing is slow when \(|\rho|\) is near \(1\).
Using the output
- Burn-in: discard the early part of the path, before the chain has forgotten its start.
- Thinning: keep every \(k\)th draw. It reduces autocorrelation but throws away information; do it only to save storage.
- Gelman-Rubin: run \(m\) chains from overdispersed starts, keep the second halves (length \(n\) each), and compare the between-chain variance \(B=\frac{n}{m-1}\sum_j(\bar\psi_{\cdot j}-\bar\psi_{\cdot\cdot})^2\) with the mean within-chain variance \(W\). Then
\[ \hat R=\sqrt{\frac{n-1}{n}+\frac1n\frac BW}, \] and values near \(1\) are required (not sufficient) for convergence.
Slow mixing is repaired by blocking (update strongly correlated components jointly from their joint conditional) or by reparameterising to near-independence. In a regression \(y_i=\alpha+\beta x_i+\varepsilon_i\), centring the covariate at \(\bar x\) makes the intercept and slope nearly uncorrelated a posteriori, and a crawling Gibbs sampler then mixes at once. Linear-model posteriors themselves are on the Regression page.
Metropolis-Hastings
Gibbs needs full conditionals that can be sampled. Metropolis-Hastings needs only the ability to evaluate the target up to a constant.
Given the current \(\theta_t\), propose \(\theta^*\sim q(\theta_t,\cdot)\) and set \(\theta_{t+1}=\theta^*\) with probability \[ \alpha(\theta_t,\theta^*)=\min\Big\{1,\ \frac{\pi(\theta^*\mid x)\,q(\theta^*,\theta_t)}{\pi(\theta_t\mid x)\,q(\theta_t,\theta^*)}\Big\}, \] otherwise \(\theta_{t+1}=\theta_t\). The unknown \(m(x)\) cancels in the ratio, so only \(L\pi\) is needed. For a symmetric (random walk) proposal the \(q\)s cancel and \(\alpha=\min\{1,\pi(\theta^*)/\pi(\theta_t)\}\); for an independence proposal \(q(\theta,\cdot)=g\) it is \(\min\{1,w(\theta^*)/w(\theta_t)\}\) with \(w=\pi/g\).
If \(\pi(\theta)P(\theta,\theta^*)=\pi(\theta^*)P(\theta^*,\theta)\) for all \(\theta,\theta^*\), then \(\pi\) is stationary for \(P\).
Integrate over \(\theta\): \(\int\pi(\theta)P(\theta,\theta^*)\,d\theta=\pi(\theta^*)\int P(\theta^*,\theta)\,d\theta=\pi(\theta^*)\).
The M-H kernel \(P(\theta,\theta^*)=q(\theta,\theta^*)\alpha(\theta,\theta^*)\) for \(\theta^*\neq\theta\), with the rejected mass \(1-\int q\alpha\) placed on \(\theta\) itself, satisfies detailed balance with respect to \(\pi\). Hence \(\pi\) is stationary.
For \(\theta=\theta^*\) there is nothing to prove. For \(\theta\neq\theta^*\), using \(c\min\{1,a/c\}=\min\{c,a\}\) for \(c>0\), \[ \pi(\theta)q(\theta,\theta^*)\alpha(\theta,\theta^*)=\min\{\pi(\theta)q(\theta,\theta^*),\ \pi(\theta^*)q(\theta^*,\theta)\}, \] which is symmetric in \((\theta,\theta^*)\) and therefore equals \(\pi(\theta^*)q(\theta^*,\theta)\alpha(\theta^*,\theta)\).
The acceptance probability is not arbitrary: requiring detailed balance and making the moves as likely as possible forces it. If \(\pi(\theta)q(\theta,\theta^*)>\pi(\theta^*)q(\theta^*,\theta)\), the move \(\theta^*\to\theta\) is under-made, so set \(\alpha(\theta^*,\theta)=1\) and solve detailed balance for \(\alpha(\theta,\theta^*)\); that is the ratio, which is less than \(1\). Detailed balance is sufficient for invariance but not necessary; non-reversible samplers exist.
Updating component \(j\) with proposal \(q(\theta,\theta^*)=\pi(\theta_j^*\mid\theta_{-j})\), \(\theta^*_{-j}=\theta_{-j}\), gives \(\alpha\equiv1\).
Write \(\pi(\theta)=\pi(\theta_{-j})\pi(\theta_j\mid\theta_{-j})\). The M-H ratio is \[ \frac{\pi(\theta_{-j})\pi(\theta_j^*\mid\theta_{-j})\cdot\pi(\theta_j\mid\theta_{-j})}{\pi(\theta_{-j})\pi(\theta_j\mid\theta_{-j})\cdot\pi(\theta^*_j\mid\theta_{-j})}=1 . \]
So samplers can be mixed component by component: Gibbs steps where a full conditional is standard, M-H steps (Metropolis-within-Gibbs) where it is not.
Tuning a random walk. With proposal \(\theta^*=\theta_t+\varepsilon\), \(\varepsilon\sim N(0,s^2)\): a tiny \(s\) accepts nearly everything but crawls, a huge \(s\) proposes into the tails and rejects nearly everything, and both give high autocorrelation. Aim for an acceptance rate around \(0.44\) in one dimension and \(0.234\) in many. Compute \(\log\alpha\) rather than \(\alpha\), since likelihoods of many observations underflow.
Graphical Models and JAGS
A directed acyclic graph (DAG) organises a model so that each full conditional involves only a few nodes. The machine-learning view of the same objects is on the Bayesian networks and Markov networks pages.
A DAG on nodes \(V\) (data and parameters alike) encodes the joint distribution \[ \pi(V)=\prod_{v\in V}\pi\big(v\mid\mathrm{pa}(v)\big), \] where \(\mathrm{pa}(v)\) are the parents of \(v\). Equivalently (the local Markov property) each node is conditionally independent of its non-descendants given its parents.
The full conditional of a node \(v\) is proportional to the product of the factors that contain \(v\): \[ \pi(v\mid V\setminus v)\propto\pi\big(v\mid\mathrm{pa}(v)\big)\prod_{w:\,v\in\mathrm{pa}(w)}\pi\big(w\mid\mathrm{pa}(w)\big). \] So it depends only on the Markov blanket of \(v\): its parents, its children and its children’s other parents. These are the neighbours of \(v\) in the conditional independence graph obtained by joining parents that share a child and dropping the arrows (moralising).
In \(\pi(v\mid V\setminus v)\propto\pi(V)\) every factor not containing \(v\) is a constant.
Lengths \(y_i\) and ages \(x_i\) of 27 dugongs follow \(y_i\sim N(\mu_i,\sigma^2)\) with the saturating growth curve \(\mu_i=\alpha-\beta\gamma^{x_i}\), \(0<\gamma<1\). The DAG has \(\alpha,\beta,\gamma,x_i\to\mu_i\to y_i\leftarrow\tau\) with \(\tau=1/\sigma^2\), the \(i\) nodes inside a plate. With vague normal priors on \(\alpha,\beta\), \(\gamma\sim\mathrm{Unif}(0.5,1)\) and a Gamma prior on \(\tau\), the full conditionals of \(\alpha\), \(\beta\) and \(\tau\) are normal, normal and Gamma, so Gibbs steps; the conditional of \(\gamma\) is \(\propto\pi(\gamma)\prod_i\exp\{-\tau(y_i-\alpha+\beta\gamma^{x_i})^2/2\}\), which is not a standard density, so \(\gamma\) gets an M-H step. Here \(\alpha\) is the adult length, \(\alpha-\beta\) the length at birth, and \(\gamma\) the rate at which the gap closes.
JAGS automates exactly this: the model file has one line per DAG node (a stochastic node with ~, a deterministic one with <-), and JAGS picks a sampler for each full conditional. The workflow is data, model, initial values for several chains, update for burn-in, coda.samples to record, then trace plots, autocorrelation plots and \(\hat R\).
Hypothesis Testing and Bayes Factors
A Bayesian compares hypotheses, or models, by their posterior probabilities. There is no privileged null and no significance level, and accepting a hypothesis means only that it is the more probable given these data and these models.
For models \(m=0,1,\dots\) with prior probabilities \(\pi(m)\), parameters \(\theta_m\), likelihoods \(L_m\) and priors \(\pi_m(\theta_m)\), the marginal likelihood of model \(m\) is \[ m_m(x)=\int L_m(x\mid\theta_m)\,\pi_m(\theta_m)\,d\theta_m , \] with every constant kept. The Bayes factor is \(B_{ij}=m_i(x)/m_j(x)\), and \[ \frac{\pi(i\mid x)}{\pi(j\mid x)}=B_{ij}\times\frac{\pi(i)}{\pi(j)},\qquad \pi(i\mid x)=\Big[1+\sum_{j\ne i}\frac{\pi(j)}{\pi(i)}B_{ji}\Big]^{-1}. \] Posterior odds are prior odds times the Bayes factor. Jeffreys’ criterion accepts \(H_0\) when its posterior odds exceed one.
Bayes’ theorem over the discrete model index gives \(\pi(i\mid x)=m_i(x)\pi(i)/\sum_jm_j(x)\pi(j)\). Dividing two such expressions cancels the denominator; dividing numerator and denominator of one of them by \(m_i(x)\pi(i)\) gives the second formula.
Let \(X_1,\dots,X_n\) be iid \(N(\theta,1)\), \(H_0:\theta=0\), \(H_1:\theta=1\), equal prior probabilities. With \(\bar X\sim N(\theta,1/n)\), \[ B_{01}=\frac{\exp(-n\bar x^2/2)}{\exp(-n(\bar x-1)^2/2)}=\exp\Big\{\frac n2(1-2\bar x)\Big\}, \] so \(H_0\) is accepted exactly when \(\bar x<\tfrac12\), the midpoint, whatever \(n\) is; the evidence for the winner grows exponentially in \(n\).
A composite hypothesis gets its own prior, and model selection is then the same calculation as estimation one level up: embed the models in an encompassing model with parameter \((m,\theta_m)\), and the joint posterior is a mixture over models with weights \(\pi(m\mid x)\). One may select the most probable model, or average over models, weighting each model’s predictions by \(\pi(m\mid x)\).
Suppose model 0 is model 1 with \(\phi=\phi_0\): \(L_0(x\mid\theta)=L_1(x\mid\theta,\phi_0)\), and the priors agree, \(\pi_0(\theta)=\pi_1(\theta\mid\phi_0)\). Then \[ B_{01}=\frac{\pi_1(\phi_0\mid x)}{\pi_1(\phi_0)}, \] the ratio of posterior to prior density of \(\phi\) at \(\phi_0\) under the larger model (the Savage-Dickey ratio).
Let \(L_1(x\mid\phi)=\int L_1(x\mid\theta,\phi)\pi_1(\theta\mid\phi)\,d\theta\). Then \(m_1(x)=\int L_1(x\mid\phi)\pi_1(\phi)\,d\phi\), so Bayes’ theorem for \(\phi\) in model 1 reads \(\pi_1(\phi_0\mid x)=L_1(x\mid\phi_0)\pi_1(\phi_0)/m_1(x)\). By nesting, \(m_0(x)=\int L_1(x\mid\theta,\phi_0)\pi_1(\theta\mid\phi_0)\,d\theta=L_1(x\mid\phi_0)\). Divide.
In the nested setting, give \(\phi\) the prior \(\mathrm{Unif}(-c,c)\) and suppose \(\int L_1(x\mid\phi)\,d\phi<\infty\). Then \(m_1(x)=\frac1{2c}\int_{-c}^cL_1(x\mid\phi)\,d\phi\to0\) as \(c\to\infty\), while \(m_0(x)\) does not involve \(c\); so \(B_{01}\to\infty\) whatever the data.
This is the Lindley (or Bartlett) paradox. In the Savage-Dickey form: once \(c\) is large the posterior density at \(\phi_0\) stops changing, but the prior density \(1/(2c)\) there goes to zero. The same holds for any scale family \(\pi_1(\phi/c)/c\). So improper priors, harmless for estimation as long as the posterior is proper, are forbidden on parameters not shared by both models: their arbitrary constant multiplies the Bayes factor. A parameter common to both models with a common prior is the exception, since its constant cancels.
- Prior Monte Carlo: \(\hat m=\frac1N\sum_iL(x\mid\theta^{(i)})\), \(\theta^{(i)}\sim\pi\). Unbiased, but most prior draws have negligible likelihood, so the variance is huge.
- Importance sampling: \(\hat m=\frac1N\sum_iL(x\mid\theta^{(i)})\pi(\theta^{(i)})/g(\theta^{(i)})\), \(\theta^{(i)}\sim g\), best when \(g\) is close to the posterior.
- Harmonic mean: \(1/\hat m=\frac1N\sum_i1/L(x\mid\theta^{(i)})\) over posterior draws, from \(E_{\pi(\theta\mid x)}[1/L]=\int\pi(\theta)/m(x)\,d\theta=1/m(x)\). Uses the MCMC output already in hand, but often has infinite variance.
A Bayes factor estimated as a ratio \(\hat m_0/\hat m_1\) of unbiased estimates is biased upwards: by Jensen, \(E[1/\hat m_1]\ge1/E[\hat m_1]\), and \(\hat m_0\), \(\hat m_1\) are computed independently.
Write an update rule as \(\pi(\theta)\mapsto\psi(\ell,\pi)\) for a loss \(\ell(x,\theta)\). The rule is coherent if updating on \(x_1\) then \(x_2\) gives the same result as on \(x_2\) then \(x_1\), and as on both at once. The only coherent rules (Bissiri, Holmes and Walker, 2016) are \(\psi(\ell,\pi)\propto e^{\ell}\pi\); with \(\ell=\log L\) this is Bayes’ theorem. Consequently, for conditionally independent blocks, \[ B_{ij}(x_1,x_2)=B_{ij}(x_2\mid x_1)\,B_{ij}(x_1), \] where the partial Bayes factor \(B_{ij}(x_2\mid x_1)\) is computed with each model’s prior first updated by \(x_1\).
The factorisation: \(m_i(x_1,x_2)=m_i(x_2\mid x_1)\,m_i(x_1)\) for each model, where \(m_i(x_2\mid x_1)=\int L_i(x_2\mid\theta)\pi_i(\theta\mid x_1)\,d\theta\) is the marginal likelihood of \(x_2\) under the updated prior. Take the ratio over \(i\) and \(j\).
The partial Bayes factor is the practical repair for the paradox: an arbitrary constant \(c\) in \(\pi_1(\phi)\) multiplies both \(B_{01}(x_1,x_2)\) and \(B_{01}(x_1)\), so it cancels in their ratio. A training sample \(x_1\) turns the vague prior into a proper, informative one, and the test is made on \(x_2\).
The Posterior as an Optimal Update
Why Bayes’ theorem, rather than some other way of combining prior and data? A modern answer: it is the solution of an optimisation problem, the same one that makes the MLE natural, with a regulariser that keeps the answer close to the prior.
For densities \(f,g\), \(\mathrm{KL}(f\Vert g)=E_f\big[\log\frac{f(X)}{g(X)}\big]=\int f\log\frac fg\). It is not symmetric, and it is \(+\infty\) if \(f>0\) on a set where \(g=0\).
\(\mathrm{KL}(f\Vert g)\ge0\), with equality if and only if \(f=g\) almost everywhere.
By Jensen for the concave \(\log\), \(-\mathrm{KL}(f\Vert g)=E_f[\log(g/f)]\le\log E_f[g/f]=\log\int_{f>0}g\le\log1=0\). Equality in Jensen forces \(g/f\) to be constant on \(\{f>0\}\), and the constant must be \(1\).
If the data are iid from \(p_{\mathrm{true}}\), then \(\mathrm{KL}(p_{\mathrm{true}}\Vert p(\cdot\mid\theta))=E[\log p_{\mathrm{true}}(X)]-E[\log p(X\mid\theta)]\). The first term is free of \(\theta\), and the second is estimated by \(\frac1n\sum_i\log p(x_i\mid\theta)\). So maximising the likelihood minimises the estimated divergence from the truth, whether or not the model is correct; with a correct model the minimiser is \(\theta_0\).
Over all distributions \(q\) on \(\Theta\), \[ \arg\min_q\Big\{E_{q}\Big[-\sum_{i=1}^n\log p(x_i\mid\theta)\Big]+\mathrm{KL}(q\Vert\pi)\Big\}=\pi(\theta\mid x)=\frac{\prod_ip(x_i\mid\theta)\,\pi(\theta)}{Z}. \]
Let \(S(\theta)=\sum_i\log p(x_i\mid\theta)\) and \(Z=\int e^{S}\pi\,d\theta\). The objective is \[ \int q\Big[-S+\log\frac q\pi\Big]d\theta=\int q\log\frac{q}{\pi e^{S}}\,d\theta=\int q\log\frac{q}{\pi e^{S}/Z}\,d\theta-\log Z=\mathrm{KL}\Big(q\,\Big\Vert\,\frac{\pi e^S}{Z}\Big)-\log Z . \] The last term does not involve \(q\), and by Gibbs’ inequality the divergence is minimised, at zero, exactly when \(q=\pi e^S/Z\), the posterior.
The first term rewards fit and on its own would put all mass on the MLE. The second is a regulariser: it is infinite if \(q\) puts mass where \(\pi\) has none, and grows without bound as \(q\) shrinks to a point mass, so the optimum stays spread out. Replacing \(-\log p\) by another loss, or restricting \(q\) to a tractable family, gives generalised Bayesian inference and variational inference respectively.
Hierarchical Models
For data \(y_{ij}\) (observation \(i\) in group \(j\)), \[ y_{ij}\mid\theta_j\sim f(y\mid\theta_j),\qquad \theta_j\mid\phi\sim\pi(\theta\mid\phi),\qquad \phi\sim\pi(\phi). \] The \(\phi\) are hyperparameters with a hyper-prior. The posterior is proportional to the product of all three levels, \(\pi(\theta,\phi\mid y)\propto\prod_{i,j}f(y_{ij}\mid\theta_j)\prod_j\pi(\theta_j\mid\phi)\,\pi(\phi)\).
Hospital \(i\) has \(r_i\) deaths in \(n_i\) operations, \(r_i\sim\mathrm{Bin}(n_i,p_i)\); in total 208 deaths in 2814 operations, and hospital 1 is the \(0/47\) hospital from the first section. Three structures:
- Pooled: \(p_i=p\) for all \(i\); with a \(\mathrm{Beta}(a,b)\) prior, \(p\mid r\sim\mathrm{Beta}(\sum r_i+a,\ \sum(n_i-r_i)+b)\). Ignores real differences between hospitals.
- Separate: independent \(p_i\sim\mathrm{Beta}(a,b)\), so \(p_i\mid r\sim\mathrm{Beta}(r_i+a,\ n_i-r_i+b)\) independently. Hospital 1’s estimate then rests on 47 operations alone.
- Hierarchical: \(p_i\mid a,b\sim\mathrm{Beta}(a,b)\) with \((a,b)\) unknown and given a hyper-prior. Integrating out each \(p_i\) with the Beta integral,
\[ \pi(a,b\mid r)\propto\pi(a,b)\prod_{i=1}^{12}\frac{B(a+r_i,\ b+n_i-r_i)}{B(a,b)},\qquad p_i\mid a,b,r\sim\mathrm{Beta}(a+r_i,\ b+n_i-r_i). \] The data choose how similar the hospitals are. Small hospitals are pulled most strongly towards the common rate of about 7%, and hospital 1 no longer looks implausibly safe.
The justification for putting a common prior on the \(\theta_j\) is a judgement about observables rather than about parameters.
Random variables \(X_1,\dots,X_n\) are exchangeable if \(f(x_1,\dots,x_n)=f(x_{\sigma(1)},\dots,x_{\sigma(n)})\) for every permutation \(\sigma\): the labels carry no information. An infinite sequence is exchangeable if every finite stretch is.
Iid implies exchangeable, but not conversely. Drawing without replacement from an urn with one red and two white balls gives exchangeable indicators (each ordering of one red has probability \(\tfrac13\)) that are dependent, with \(\mathrm{Cov}(X_1,X_2)=-\tfrac19\). A Markov chain, even a stationary one, is generally not exchangeable.
If \(X_1,X_2,\dots\) is an infinitely exchangeable sequence, then for every \(n\) \[ f(x_1,\dots,x_n)=\int\prod_{i=1}^nf(x_i\mid\theta)\,\pi(\theta)\,d\theta \] for some parameter \(\theta\), sampling model \(f(\cdot\mid\theta)\) and distribution \(\pi\).
So exchangeable quantities behave as iid draws given an unknown parameter that itself has a prior: the hierarchical model, derived from a judgement of similarity. Infinite exchangeability is essential. The three-ball urn cannot be represented this way, since any such mixture has \(P(X_1=X_2=1)=E[\theta^2]\ge(E\theta)^2=\tfrac19\), whereas the urn gives \(0\).
If \(y_i\mid\theta_i\sim N(\theta_i,s_i^2)\) and \(\theta_i\sim N(\mu,\sigma^2)\) with \(\mu,\sigma^2,s_i^2\) known, then \[ \theta_i\mid y_i\sim N\big(w_i\mu+(1-w_i)y_i,\ (1-w_i)s_i^2\big),\qquad w_i=\frac{s_i^2}{s_i^2+\sigma^2}. \]
This is the known-variance normal update with prior precision \(1/\sigma^2\), one observation of precision \(1/s_i^2\): the posterior precision is \(\frac1{s_i^2}+\frac1{\sigma^2}=\frac{s_i^2+\sigma^2}{s_i^2\sigma^2}\), so the variance is \(\frac{\sigma^2}{s_i^2+\sigma^2}s_i^2=(1-w_i)s_i^2\), and the mean is the precision-weighted average \(\frac{y_i/s_i^2+\mu/\sigma^2}{1/s_i^2+1/\sigma^2}=w_i\mu+(1-w_i)y_i\).
Each estimate is pulled towards the common mean, and its interval narrows, by the same factor \(w_i\). Noisy units (large \(s_i^2\)) shrink most; when units vary little (small \(\sigma^2\)) everything shrinks to the common mean.
Eight schools report estimated coaching effects \(\bar y_{\cdot j}\) with known standard errors \(\sigma_j\). Model \(\bar y_{\cdot j}\sim N(\theta_j,\sigma_j^2)\), \(\theta_j\sim N(\mu,\sigma^2)\), \(\pi(\mu)\propto1\). Given the between-school variance \(\sigma^2\), integrating out \(\theta_j\) gives \(\bar y_{\cdot j}\mid\mu,\sigma^2\sim N(\mu,\sigma_j^2+\sigma^2)\), so \(\mu\mid\sigma^2,y\) is normal around the precision-weighted mean \(\hat\mu=\sum_j\frac{\bar y_{\cdot j}}{\sigma_j^2+\sigma^2}\big/\sum_j\frac{1}{\sigma_j^2+\sigma^2}\), and \[ E[\theta_j\mid\sigma^2,y]=w_j\hat\mu+(1-w_j)\bar y_{\cdot j},\qquad w_j=\frac{\sigma_j^2}{\sigma_j^2+\sigma^2}. \] As \(\sigma\to0\) every school gets the pooled estimate; as \(\sigma\to\infty\) each keeps its own. The hierarchical model puts a posterior on \(\sigma\) and lands in between: the schools borrow strength from one another.
Mixture Models
An \(M\)-component mixture has density \[ f(x\mid\theta,\lambda)=\sum_{m=1}^M\lambda_mf_m(x\mid\theta_m),\qquad\lambda_m\ge0,\quad\sum_m\lambda_m=1, \] describing a population of unlabelled sub-populations.
The likelihood \(\prod_i\sum_m\lambda_mf_m(x_i\mid\theta_m)\) has \(M^n\) terms when expanded, and no conjugate prior. The trick is to add the missing labels as parameters.
Let \(z_i=(z_{i1},\dots,z_{iM})\) be one-hot with \(z_i\mid\lambda\sim\mathrm{Multinomial}(1;\lambda)\), and \(x_i\mid z_{im}=1\sim f_m(\cdot\mid\theta_m)\). Then \[ f(x,z\mid\theta,\lambda)=\prod_{i=1}^n\prod_{m=1}^M\big[\lambda_mf_m(x_i\mid\theta_m)\big]^{z_{im}}, \] and summing over each \(z_i\) recovers the mixture likelihood.
With a \(\mathrm{Dirichlet}(\alpha_1,\dots,\alpha_M)\) prior on \(\lambda\) and independent priors on the \(\theta_m\), the full conditionals are
- \(\theta_m\mid\cdots\): the one-population posterior using only the observations currently allocated to component \(m\);
- \(\lambda\mid\cdots\sim\mathrm{Dirichlet}(\alpha_1+n_1,\dots,\alpha_M+n_M)\), where \(n_m=\sum_iz_{im}\);
- \(z_i\mid\cdots\sim\mathrm{Multinomial}(1;\tilde p_{i1},\dots,\tilde p_{iM})\) independently, with \(\tilde p_{im}\propto\lambda_mf_m(x_i\mid\theta_m)\).
In the joint \(\pi(\theta)\pi(\lambda)f(x,z\mid\theta,\lambda)\), the factors containing \(\theta_m\) are \(\pi(\theta_m)\prod_{i:z_{im}=1}f_m(x_i\mid\theta_m)\); those containing \(\lambda\) are \(\prod_m\lambda_m^{\alpha_m-1}\cdot\prod_m\lambda_m^{n_m}\); and those containing \(z_i\) are \(\prod_m[\lambda_mf_m(x_i\mid\theta_m)]^{z_{im}}\), which normalised over the \(M\) one-hot values is the stated multinomial.
If a component is empty (\(n_m=0\)) its full conditional is its prior, so the priors on the \(\theta_m\) must be proper. The posterior probability that observation \(i\) belongs to component \(m\) is estimated by averaging \(\tilde p_{im}\) (or \(z_{im}\)) over the draws; a new observation is predicted by composition: draw \((\theta,\lambda)\) from the posterior, then a component \(m^*\sim\lambda\), then \(y^*\sim f_{m^*}(\cdot\mid\theta_{m^*})\). For two equal-variance normal components the plug-in classification probability is a logistic curve in \(x\). The EM algorithm for the same model is on the Gaussian mixture models page.
Permuting the component labels of \((\lambda_m,\theta_m,z_{im})\) leaves the likelihood unchanged, so the mixture is not identifiable: \(M!\) parameter values give the same model. Under a prior symmetric in the components, every \(\theta_m\) has the same marginal posterior.
A permutation \(s\) only reorders the finite sum \(\sum_m\lambda_{s(m)}f(x\mid\theta_{s(m)})=\sum_k\lambda_kf(x\mid\theta_k)\). If the prior is also invariant, so is the posterior, and the swap \(1\leftrightarrow2\) maps \(\{\theta_1\in A\}\) to \(\{\theta_2\in A\}\) with equal probability.
A sampler that mixes well therefore produces marginal posteriors that are all the same useless average of the components. The fix is an identifying constraint (order the means, or the weights) or an informative asymmetric prior, applied before the components are summarised.
For yearly counts \(y_1,\dots,y_n\), let \(y_i\sim\mathrm{Pois}(\theta_1)\) for \(i\le\tau\) and \(\mathrm{Pois}(\theta_2)\) after, with \(\theta_1,\theta_2\sim\mathrm{Gamma}(a,b)\) and \(\tau\) discrete uniform. This is a mixture whose indicators are all fixed by the single parameter \(\tau\). With \(S_k=\sum_{i\le k}y_i\), the Gibbs updates are \[ \theta_1\mid\tau,y\sim\mathrm{Gamma}(a+S_\tau,\ b+\tau),\qquad \theta_2\mid\tau,y\sim\mathrm{Gamma}(a+S_n-S_\tau,\ b+n-\tau), \] and \(\tau\) is drawn from its discrete full conditional \(\pi(\tau\mid\theta,y)\propto\theta_1^{S_\tau}e^{-\tau\theta_1}\,\theta_2^{S_n-S_\tau}e^{-(n-\tau)\theta_2}\), computed for every \(\tau\) and normalised. Whether there is a change at all is a Bayes factor between this model and a constant rate, which needs proper priors.
Mixtures can carry a scientific hypothesis inside a hierarchical model. In a reaction-time study (from Gelman et al., Bayesian Data Analysis, chapter 18), 17 subjects, six of them schizophrenic, each give 30 log reaction times. The theory says schizophrenics have a general motor slowing (a shift in every trial) and an attentional deficit on some trials (a delay \(\tau\) on a random subset). So each subject has a random mean \(\alpha_i\sim N(\mu+\beta S_i,\sigma_\alpha^2)\), with \(S_i\) the diagnosis, and trial \(j\) is delayed by \(\tau\) when an indicator \(z_{ij}\sim\mathrm{Bernoulli}(\lambda S_i)\) is one. A posterior predictive check (simulate replicated data from the posterior and compare statistics aimed at the doubtful assumption, here the smallest and largest within-person standard deviations, with their observed values) shows this is too rigid, and the expanded model lets only a proportion \(\omega\) of schizophrenics have delays, with extra variance on the delayed trials. The number of components can also be chosen by posterior model probabilities, as in the Bayes factor section.
Where It Leads
- Better samplers: Hamiltonian Monte Carlo, adaptive MCMC, reversible-jump MCMC for models of varying dimension, and sequential Monte Carlo (particle methods). The Approximate inference page covers sampling in graphical models.
- Approximations: variational Bayes restricts \(q\) in the optimal-update problem to a tractable family; approximate Bayesian computation replaces a likelihood that cannot be evaluated by simulation from the model.
- Generalised Bayes: replace \(-\log p\) in the optimal-update objective by a loss tailored to the question, keeping coherence.
- Foundations: the posterior is a conditional distribution, and its rigorous definition (regular conditional probabilities, the Radon-Nikodym derivative) is in Measure Theory.
Backlinks (1)
1. Statistics /wiki/mathematics/statistics/
This page pairs well with Probability.
Two child pages go further:
- Bayesian Inference: priors and posteriors, decisions, Monte Carlo and MCMC, Bayes factors, hierarchical and mixture models.
- Multivariate Analysis: the multivariate normal and Wishart distributions, Hotelling’s \(T^2\), tests on mean vectors and covariances, PCA, MANOVA, factor analysis, classification and copulas.
Table of Distributions
| Distribution | mass/density function | \[S_X\] | \[\mathbb{E}(X)\] | \[\mathrm{Var}(X)\] | \[\phi_X(s)\] |
|---|---|---|---|---|---|
| Bernoulli \[\mathrm{Bern}(\pi)\] | \[P(X=1)=\pi\\ P(X=0)=1-\pi\] | \[\{0,1\}\] | \[\pi\] | \[\pi(1-\pi)\] | \[(1-\pi)+\pi e^{s}\] |
| Binomial \[\mathrm{Bin}(n,\pi)\] | \[p_X(x)=\binom{n}{x}\pi^{x}(1-\pi)^{n-x}\] | \[\{0,1,\dots,n\}\] | \[n\pi\] | \[n\pi(1-\pi)\] | \[(1-\pi+\pi e^{s})^{n}\] |
| Geometric \[\mathrm{Geo}(\pi)\] | \[p_X(x)=\pi(1-\pi)^{x-1}\] | \[\{1,2,\dots\}\] | \[\pi^{-1}\] | \[(1-\pi)\pi^{-2}\] | \[\frac{\pi}{e^{-s}-1+\pi}\] |
| Poisson \[\mathcal{P}(\lambda)\] | \[p_X(x)=e^{-\lambda}\lambda^{x}/x!\] | \[\{0,1,\dots\}\] | \[\lambda\] | \[\lambda\] | \[\exp\{\lambda(e^{s}-1)\}\] |
| Uniform \[U[\alpha,\beta]\] | \[f_X(x)=(\beta-\alpha)^{-1}\] | \[[\alpha,\beta]\] | \[\frac{1}{2}(\alpha+\beta)\] | \[\frac{1}{12}(\beta-\alpha)^2\] | \[\frac{e^{\beta s}-e^{\alpha s}}{s(\beta-\alpha)}\] |
| Exponential \[\mathrm{Exp}(\lambda)\] | \[f_X(x)=\lambda e^{-\lambda x}\] | \[[0,\infty)\] | \[\lambda^{-1}\] | \[\lambda^{-2}\] | \[\frac{\lambda}{\lambda-s}\] |
| Gaussian \[\mathcal{N}(\mu,\sigma^{2})\] | \[f_X(x)=\frac{1}{\sqrt{2\pi\sigma^{2}}}\exp\left\{-\frac{(x-\mu)^2}{2\sigma^{2}}\right\}\] | \[\mathbb{R}\] | \[\mu\] | \[\sigma^{2}\] | \[e^{\mu s+\frac{1}{2}\sigma^{2}s^{2}}\] |
| Gamma \[\Gamma(\alpha,\lambda)\] | \[f_X(x)=\frac{1}{\Gamma(\alpha)}\lambda^{\alpha}x^{\alpha-1}e^{-\lambda x}\] | \[[0,\infty)\] | \[\alpha\lambda^{-1}\] | \[\alpha\lambda^{-2}\] | \[\left(\frac{\lambda}{\lambda-s}\right)^{\alpha}\] |
Statistical Inference
Let \(X=(X_1,\ldots,X_n)\) be i.i.d. from a parametric family \(\{F_\theta:\theta\in\Theta\subset\mathbb{R}^p\}\). The parameter \(\theta\) is unknown; inference uses the randomness of \(X\) to learn about \(\theta\).
Comments