Chapter 05: Bayesian Statistics

This note is a high-fidelity Markdown migration of the Bayesian Statistics chapter from the LaTeX source.

Parent map: index Prerequisites: probability-and-mathstats, linear-regression, maximum-likelihood-and-machine-learning

Concept map

flowchart TD
  A[Bayes Theorem] --> B[Priors]
  B --> C[Conjugate Updating]
  A --> D[Posterior Predictive]
  A --> E[Model Selection]
  E --> F[Marginal Likelihood]
  C --> G[Hierarchical Models]
  G --> H[Empirical Bayes]
  A --> I[Computation]
  I --> J[MCMC]
  J --> K[Gibbs]
  J --> L[Metropolis-Hastings]
  I --> M[EM Algorithm]
  A --> N[Graphical Models]

Bayesian Statistics

Setup

Notation: per the Murphy textbook, some statements use notation $\mathcal{D}:={\left\{({\mathbf{x}}_i,y_i)\right\}}_{i=1}^N$ as shorthand for data.

$$f(\mathbf{\theta }\mid \mathbf{X})=\frac{f(\mathbf{X}\mid \mathbf{\theta })f(\mathbf{\theta })}{\int f(\mathbf{X}\mid \mathbf{\theta })f(\mathbf{\theta })d\mathbf{\theta }}\propto \mathcal{L}(\mathbf{\theta })f(\mathbf{\theta })$$

1. Use Bayes Rule to come up with a posterior probability of some hypothesis $H_u$ given event $E$. Your prior is $P(H_u)$.

$$\begin{array}{r}P(H_u|E)=\frac{P(E|H_u)P(H_u)}{P(E|H_u)P(H_u)+P(E|H_u^c)P(H_u^c)}\end{array}$$

Call the posterior probability $P(H_u^{\prime })$.

2. Given second event $E^{\prime }$, use posterior probability from step 1 as your prior in the second update step.

$$\begin{array}{r}P(H_u|E^{\prime })=\frac{P(E^{\prime }|H_u)P(H_u^{\prime })}{P(E^{\prime }|H_u)P(H_u^{\prime })+P(E^{\prime }|H_u^c)P((H_u^{\prime }{)}^c)}\end{array}$$

A sequence of random variables $y_1,\dots ,y_n$ is finitely exchangeable if their joint density remains the same under any re-ordering or re-labeling of the indices of the data.

$$p(y_1,\dots ,y_n)=p(y_{z(1)},\dots ,y_{z(n)}).$$

Exchangeability justifies use of the prior: If the data are exchangeable, then there is a parameter $\theta$ that drives the stochastic model generating the data and there exists a density over $\theta$ that does not depend on the data itself. The data are conditionally i.i.d., given the prior $\theta$. Independence vs. Exchangeability: Independence is a stronger condition than exchangeability (it is a special case of exchangeability). Exchangeability only requires that the marginal distribution of each random variable is the same, i.e. $p(y_1)=p(y_2)$. Independence requires that $p(y_1|y_2)=p(y_1)$. As a result, you can have exchangeability in situations where you do not have independence, most notably sampling without replacement. If the marginal probabilities are unknown, then we only have exchangeability (not independence) even if the samples are drawn with replacement, due to the possibility that there is only one unit with a particular value of $y$.

With the full posterior, one can compute Posterior Mean, median, and mode (the latter is sometimes called the Maximum A Posteriori estimate).

One can also compute a Highest Posterior Density Region $R(\theta )$, which is a region such that the parameter lies in the region with probability $1-\alpha$.

$$1-\alpha =\mathbf{Pr}\left(\theta \in R(\theta )|y\right)={\int }_{R(\theta )}p(\theta |y)d\theta$$

Consider out-of-sample prediction for a single observation $\tilde{y}$. The posterior predictive density is

$$\begin{array}{r}p(\tilde{y}|y_1,\dots ,y_n)={\int }_{-\infty }^{\infty }p(\tilde{y}|\theta ,y_1,\dots ,y_n)p(\theta |y_1,\dots ,y_n)d\theta \end{array}$$

Because $\tilde{y}$ is independent of $y$ conditional on $\theta$ (exchangeability), we can simplify this as

$$\begin{array}{rl}p(\tilde{y}|y_1,\dots ,y_n)& ={\int }_{-\infty }^{\infty }p(\tilde{y}|\theta ,y_1,\dots ,y_n)p(\theta |y_1,\dots ,y_n)d\theta \\ & ={\int }_{-\infty }^{\infty }p(\tilde{y}|\theta )p(\theta |y_1,\dots ,y_n)d\theta \end{array}$$

This is just the data density for $y$ multiplied by the posterior density for $\theta$.

Consider $\tilde{y}\sim \text{Bernoulli}(\theta )$. The posterior predictive density is

$$\begin{array}{rl}p(\tilde{y}|y)& ={\int }_0^{1}p(\tilde{y}|\theta )p(\theta |y)d\theta \\ & ={\int }_0^1{\theta }^{\tilde{y}}(1-\theta {)}^{1-\tilde{y}}p(\theta |y)d\theta \end{array}$$

So if we want to know the posterior predictive probability $p(\tilde{y}=1|\theta )$, we can compute it as

$$\begin{array}{rl}p(\tilde{y}=1|\theta )& ={\int }_0^1\theta p(\theta |y)d\theta \\ & =E[\theta |y]\end{array}$$

which is the posterior mean.

An uninformative prior on $\theta$ produces a posterior density that is proportional to the likelihood (differing only by the constant of proportionality). This implies that the mode of the posterior density is the $\theta$ that maximizes the likelihood function. An informative prior on $\theta$ yields a posterior mean that is a precision-weighted average of the prior mean and the MLE.

Stan dev team recommendations: https://github.com/stan-dev/stan/wiki/Prior-Choice-Recommendations

$$\theta |y{\sim }^a\mathcal{N}\left(\hat{\theta },\mathcal{I}(\hat{\theta }{)}^{-1}\right)$$

As $N\to \infty$, the likelihood component of the posterior becomes dominant and as a result frequentist and bayesian inferences will be based on the same limiting multivariate normal distribution.

To choose between Bayesian models, we compute the posterior over models

$$p(m|\mathcal{D})=\frac{p(\mathcal{D}|m)p(m)}{{\sum }_{m\in \mathcal{M}}p(m,\mathcal{D})}$$

which allows us to pick the MAP model $\hat{m}=\mathrm{argmax}p(m|\mathcal{D})$. If we use a uniform prior over models $p(m)\propto 1$, this amounts to picking the model wich maximises

$$p(\mathcal{D}|m)=\int p(\mathcal{D}|\mathbf{\theta })p(\mathbf{\theta }|m)d\mathbf{\theta }$$

which is called the marginal likelihood / integrated likelihood / evidence for model $m$.

Conjugate Priors and Updating

Priors of the form $f(\theta )\propto cc>0$ are improper because $\int f(\theta )d\theta =\infty$. Improper priors generally not a problem as long as resulting posterior is well defined.

Suppose $X\sim \mathrm{Bernoulli}\left(p\right)$, and we choose prior $f(p)$ = 1. Define transformation $\psi =\log (p/(1-p))$. Resulting distribution of $\psi$ is $f_{\psi }(\psi )=\frac{e^{\psi }}{(1+e^{\psi }{)}^2}$

Jeffreys’ Prior is a method of constructing invariant priors.

$f(\theta )\propto I(\theta {)}^{1/2}$. For multiparameter model, $f(\theta )\propto {\left|I(\theta )\right|}^{1/2}$

Given $\gamma =h(\theta )$, $\frac{\partial \mathcal{L}}{\partial \gamma }=\frac{\partial \mathcal{L}}{\partial \theta }=\frac{\partial \theta }{\partial \gamma }$ and

$$\frac{{\partial }^2\mathcal{L}}{\partial {\gamma }^2}=\frac{{\partial }^2\mathcal{L}}{\partial {\theta }^2}{\left(\frac{\partial \theta }{\partial \gamma }\right)}^2+\frac{\partial \mathcal{L}}{\partial \theta }\frac{{\partial }^2\theta }{\partial {\gamma }^2}$$

Taking expectations wrt sample density sends second piece to zero (since $\mathbb{E}\left[\frac{\partial \mathcal{L}}{\partial \theta }\right]=0$), so

$$\mathcal{I}(\gamma )=\mathcal{I}(\theta ){\left(\frac{\partial \theta }{\partial \gamma }\right)}^2⟹{\left|\mathcal{I}(\gamma )\right|}^{1/2}={\left|\mathcal{I}(\theta )\right|}^{1/2}\left|\frac{\partial \theta }{\partial \gamma }\right|$$

Analytically tractable expressions for the posterior are derived when sample and prior densities form a natural conjugate pair, defined as having the property that sample, prior, and posterior densities all lie in the same class of densities.

Exponential family is essentially the only class of densities to have natural conjugate priors.

A one parameter member of the exponential family has density for N obs that can be expressed as

$$\mathcal{L}(y|\theta )=\prod \exp \left((a(\theta ))+b(y)+c\left(\theta \right)u\left(y\right)\right)\propto \exp \left(Na(\theta )+c(\theta )\sum _{i}u(y)\right)$$

Let $X_1,\dots ,X_n\sim \mathrm{Bernoulli}\left(p\right)$, and we take prior $f(p)=1$. By Bayes thm, the posterior is of the form

$$f(p|x^n)\propto f(p){\mathcal{L}}_n(p)=p^s(1-p{)}^{n-s}=p^{\sum x_i}(1-p{)}^{n-\sum x_i}$$

Instead we take $f(p)=\mathrm{Beta}\left(\alpha ,\beta \right)$. Uniform prior is a special case with $\alpha =\beta =1$. In general, the posterior is of the form

$$\begin{array}{rl}f(p|x^n)& =\frac{\Gamma (n+2)}{\Gamma (s+1)\Gamma (n-s+1)}{p}^{\sum x_i}(1-p{)}^{n-\sum x_i}\\ p|x^n& \sim \mathrm{Beta}\left(\sum x_i+1,n-\sum x_i+1\right)\\ & \sim \mathrm{Beta}\left({\alpha }^{\prime },{\beta }^{\prime }\right)\\ & \text{rbeta(n, shape1, shape2)}\end{array}$$

Quantity	Formula
Posterior Mean	$\frac{{\alpha }^{\prime }}{{\alpha }^{\prime }+{\beta }^{\prime }}=\frac{\sum x_i+1}{n+2}$
Posterior Mode	$\frac{{\alpha }^{\prime }-1}{{\alpha }^{\prime }+{\beta }^{\prime }-2}=\frac{\sum x_i}{n}$
Posterior Variance	$\frac{{\alpha }^{\prime }{\beta }^{\prime }}{({\alpha }^{\prime }+{\beta }^{\prime }{)}^2({\alpha }^{\prime }+{\beta }^{\prime }+1)}$
Posterior Predictive Distribution	$\sim \mathrm{Beta}\text{-}\mathrm{Binomial}(n,a,b)$ with updated parameters $\mathrm{Beta}\text{-}\mathrm{Binomial}(n,\alpha +\sum x_i,\beta +n-\sum x_i)$; library(extraDistr), rbbinom(n, size, alpha, beta)

Suppose we have a proportion from a previous study ${\theta }_0$, with variance $V({\theta }_0;\alpha ,\beta )$. Then we can create a constant

$$\begin{array}{r}\gamma =\frac{{\theta }_0(1-{\theta }_0)}{V({\theta }_0;\alpha ,\beta )}-1\end{array}$$

And compute the hyper-parameters $\alpha$ and $\beta$ for our prior distribution as

$$\begin{array}{rl}\alpha & =\gamma {\theta }_0\\ \beta & =\gamma (1-{\theta }_0)\end{array}$$

Surprisingly this works! See Jackman p.55 for a worked out example.

Let $Y_1,\dots ,Y_n\sim \text{Poisson}(\lambda )$. This means that

$$\begin{array}{rl}p(\mathbf{y}|\lambda )& =\prod _{i=1}^N\frac{\exp (-\lambda ){\lambda }^{y_i}}{y_i!}\\ & \propto {\lambda }^{\sum y_i}\exp (-n\lambda )\end{array}$$

We specify a Gamma prior on $\lambda$, which has density

$$\begin{array}{rl}p(\lambda ;a,b)& =\frac{b^a}{\Gamma (a)}{\lambda }^{a-1}\exp (-b\lambda )\end{array}$$

So then the posterior for $\lambda$ is

$$\begin{array}{rl}p(\lambda |\mathbf{y})& \propto p(\lambda )p(\mathbf{y}|\lambda )\\ & \propto {\lambda }^{a-1}\exp (-b\lambda ){\lambda }^{\sum y_i}\exp (-n\lambda )\\ & ={\lambda }^{\sum y_i+a-1}\exp (-\lambda (b+n))\\ & \sim \text{Gamma}(\sum y_i+a,b+n)\\ & \sim \text{Gamma}(a^{\prime },b^{\prime })\\ & \text{rgamma(n, shape, rate)}\end{array}$$

A flat prior is $a=b=0$.

Quantity	Formula
Posterior Mean	$\frac{a^{\prime }}{b^{\prime }}=\frac{\sum y_i+a}{b+n}$
Posterior Mode	$\frac{a^{\prime }-1}{b^{\prime }}=\frac{\sum y_i+a-1}{b+n}$
Posterior Variance	$\frac{a^{\prime }}{(b^{\prime }{)}^2}=\frac{\sum y_i+a}{(b+n{)}^2}$
Posterior Predictive Distribution	$\sim \mathrm{Negative}\mathrm{Binomial}(y,\theta )$ and $\sim \mathrm{Negative}\mathrm{Binomial}\left(a,1-\frac{1}{b+1}\right)$; rnbinom(n, size, prob)

$$\begin{array}{rl}p(\theta |{\alpha }_1,\dots {\alpha }_k)& \propto \prod _{j=1}^K{\theta }_j^{{\alpha }_j-1}\\ p(y|\theta )& \propto \prod _{j=1}^K{\theta }_j^{y_j}\end{array}$$

where $y_j$ is the count of observations in category $j$. For 3 categories, the posterior is:

$$\begin{array}{rl}p({\theta }_1,{\theta }_2,1-{\theta }_1-{\theta }_2|y)& \propto {\theta }_1^{{\alpha }_1+y_1-1}{\theta }_2^{{\alpha }_2+y_2-1}(1-{\theta }_1-{\theta }_2{)}^{{\alpha }_3+y_3-1}\\ & \sim \text{Dirichlet}({\alpha }_1+y_1,{\alpha }_2+y_2,{\alpha }_3+y_3)\end{array}$$

$y\sim \mathcal{N}\left(\mu ,{\sigma }^2\right)$, where ${\sigma }^2$ is known but mean $\mu$ is not known. The joint density of $y$ is

$$\mathcal{L}(y|\theta )=\prod _{i=1}^N(2\pi {\sigma }^2{)}^{-1/2}\exp \left(-\frac{(y_i-\theta {)}^2}{2{\sigma }^2}\right)\propto \exp \left(-\frac{N}{2{\sigma }^2}(\bar{y}-\theta {)}^2\right)$$

Given a normal prior $\theta \sim \mathcal{N}\left(\mu ,{\tau }^2\right)⟹f(\theta )\propto \exp \left(-\frac{(\theta -\mu {)}^2}{2{\tau }^2}\right)$, we can write the posterior density of the form

$$f(\theta |y)\propto \exp \left(-\frac{N}{2{\sigma }^2}(\theta -\bar{y}{)}^2\right)\exp \left(-\frac{(\theta -\mu {)}^2}{2{\tau }^2}\right)\propto \exp \left(-\frac{1}{2}\left(\frac{(\theta -{\mu }_1{)}^2}{{\tau }_1^2}\right)\right)$$

where ${\mu }_1={\tau }_1^2(N\bar{y}/{\sigma }^2+\mu /{\tau }^2)$ and ${\tau }_1^2=(N/{\sigma }^2+1/{\tau }^2{)}^{-1}$. Posterior mean is a weighted sum of prior mean $\mu$ and sample mean $\bar{y}$ with weights that reflect the precision of the likelihood via $N/{\sigma }^2$ and prior precision $1/{\tau }^2$. Three cases (ref. Jackman p.80-94):

1. Variance known, mean unknown. Model:

$$\begin{array}{r}\mu \sim \mathcal{N}({\mu }_0,{\sigma }_0^2)\\ y\sim \mathcal{N}(\mu ,{\sigma }^2)\end{array}$$

Quantity	Formula
Posterior Mean	$\frac{{\mu }_0(1/{\sigma }_0^2)+\bar{y}(n/{\sigma }^2)}{(1/{\sigma }_0^2)+(n/{\sigma }^2)}=\bar{y}\frac{n{\sigma }_0^2}{n{\sigma }_0^2+{\sigma }^2}+{\mu }_0\frac{{\sigma }^2}{n{\sigma }_0^2+{\sigma }^2}$
Posterior Variance	${\left(\frac{1}{{\sigma }_0^2}+\frac{n}{{\sigma }^2}\right)}^{-1}=\frac{{\sigma }_0^2{\sigma }^2}{{\sigma }^2+n{\sigma }_0^2}$
Posterior Predictive Distribution	$\sim \mathcal{N}(\tilde{\mu },{\tilde{\sigma }}^2)$, with $\tilde{\mu }=\frac{n_0{\mu }_0+n\bar{y}}{n_0+n}$ and ${\tilde{\sigma }}^2={\sigma }^2+{\left(\frac{1}{{\sigma }_0^2}+\frac{n}{{\sigma }^2}\right)}^{-1}$

2. Variance and mean both unknown. Prior densities:

$$\begin{array}{rl}p(\mu ,{\sigma }^2)& =p(\mu |{\sigma }^2)p({\sigma }^2)\\ p(\mu |{\sigma }^2)& \sim \text{Normal}({\mu }_0,{\sigma }^2/n_0)\\ p({\sigma }^2)& \sim \text{Scaled-Invese-χ2}({\nu }_0/2,{\nu }_0{\sigma }_0^2/2)\end{array}$$

Conditional posterior densities: 

$$\begin{array}{rl}\mu |{\sigma }^2,y& \sim \mathcal{N}({\mu }_1,{\sigma }^2/n_1)\\ {\sigma }^2|y& \sim \text{Scaled-Invese-χ2}({\nu }_1/2,{\nu }_1{\sigma }_1^2/2)\end{array}$$

where

$$\begin{array}{rl}{n}_1& =n_0+n\\ {\mu }_1& =\frac{n_0{\mu }_0+n\bar{y}}{n_1}\\ {\nu }_1& ={\nu }_0+n\\ {\nu }_1{\sigma }_1^2& =v_0{\sigma }_0^2+\sum _{i=1}^N(y_i-\bar{y}{)}^2+\frac{n_{0}n}{n_0+n}({\mu }_0-\bar{y}{)}^2\end{array}$$

Marginal posterior density of $\mu$: 

$$\begin{array}{rl}p(\mu )& \sim \text{Student-T}({\mu }_1,\sqrt{{\sigma }_1^2/n_1},v_1)\\ & \sim \text{brms::rstudent\_t(n, df, mu = 0, sigma = 1)}\end{array}$$

where

$$\begin{array}{rl}{n}_1& =n_0+n\\ {\mu }_1& =\frac{n_0{\mu }_0+n\bar{y}}{n_1}\\ {\nu }_1& ={\nu }_0+n\\ {\sigma }_1^2& =S_1/{\nu }_1\\ S_1& ={\nu }_0{\sigma }_0^2+(n-1)\sum _{i=1}^N(y_i-\bar{y}{)}^2+\frac{n_{0}n}{n_1}(\bar{y}-{\mu }_0{)}^2\end{array}$$

Posterior predictive distribution for $\tilde{y}$: 

$$\begin{array}{r}p(\tilde{y}|y)\sim \text{Student-T}({\mu }_1,{\sigma }_1\sqrt{(n_1+1)/n_1},v_1)\end{array}$$

where

$$\begin{array}{rl}{n}_1& =n_0+n\\ {\mu }_1& =\frac{n_0{\mu }_0+n\bar{y}}{n_1}\\ {\nu }_1& ={\nu }_0+n\\ {\sigma }_1^2& =S_1/{\nu }_1\\ S_1& ={\nu }_0{\sigma }_0^2+(n-1)\sum _{i=1}^N(y_i-\bar{y}{)}^2+\frac{n_{0}n}{n_1}(\bar{y}-{\mu }_0{)}^2\end{array}$$

3. Improper reference prior. Prior densities:

$$\begin{array}{r}p(\mu ,{\sigma }^2)\propto 1/{\sigma }^2\end{array}$$

Posterior densities 

$$\begin{array}{rl}\mu |{\sigma }^2,y& \sim \mathcal{N}(\bar{y},{\sigma }^2/n)\\ {\sigma }^2|y& \sim \text{Scaled-Invese-χ2}(\frac{n-1}{2},\frac{{\sum }_{i=1}^N(y_i-\bar{y}{)}^2}{2})\end{array}$$

which implies

$$\begin{array}{r}\frac{\mu -\bar{y}}{\sqrt{S/((n-1)n)}}\sim t_{n-1}\end{array}$$

Posterior predictive distribution 

$$\begin{array}{r}p(\tilde{y}|y)\sim \text{Student-T}(\bar{y},s\sqrt{\frac{n+1}{n}},n-1)\end{array}$$

where

$$\begin{array}{rl}{s}^2& =\frac{1}{n-1}\sum _{i=1}^n(y-\bar{y}{)}^2\\ \bar{y}& =\frac{1}{n}\sum _{i=1}^n{y}_i\end{array}$$

Posterior can often be approximated by simulation.

• Draw ${\theta }_1,\dots ,{\theta }_B\sim p(\theta |x^n)$.
• Histogram of ${\theta }_1,\dots ,{\theta }_B$ approximates posterior density $p(\theta |x^n)$.

Methods for this: Markov-Chain Monte-Carlo, Metropolis-Hastings, Hamiltonian Monte-Carlo.

Conjugacy for Discrete Distributions

Likelihood	Conjugate prior	Posterior hyperparameters
Bernoulli$(p)$	Beta$(\alpha ,\beta )$	$\alpha +{\sum }_i{x}_i,\beta +n-{\sum }_i{x}_i$
Binomial$(p)$	Beta$(\alpha ,\beta )$	$\alpha +{\sum }_i{x}_i,\beta +{\sum }_i{N}_i-{\sum }_i{x}_i$
Negative Binomial$(p)$	Beta$(\alpha ,\beta )$	$\alpha +rn,\beta +{\sum }_i{x}_i$
Poisson$(\lambda )$	Gamma$(\alpha ,\beta )$	$\alpha +{\sum }_i{x}_i,\beta +n$
Multinomial$(p)$	Dirichlet$(\alpha )$	$\alpha +{\sum }_i{x}^{(i)}$

Conjugacy for Continuous Distributions

Likelihood	Conjugate prior	Posterior hyperparameters
Uniform$[0,\theta ]$	Pareto$(x_m,k)$	$\max \{x_{(n)},x_m\},k+n$
Exponential$(\lambda )$	Gamma$(\alpha ,\beta )$	$\alpha +n,\beta +{\sum }_i{x}_i$
Normal$(\mu ,{\sigma }_c^2)$	Normal$({\mu }_0,{\sigma }_0^2)$	$\left(\frac{{\mu }_0}{{\sigma }_0^2}+\frac{{\sum }_i{x}_i}{{\sigma }_c^2}\right)/\left(\frac{1}{{\sigma }_0^2}+\frac{n}{{\sigma }_c^2}\right),{\left(\frac{1}{{\sigma }_0^2}+\frac{n}{{\sigma }_c^2}\right)}^{-1}$
Normal$({\mu }_c,{\sigma }^2)$	Scaled Inverse Chi-square$(\nu ,{\sigma }_0^2)$	$\nu +n,\frac{\nu {\sigma }_0^2+{\sum }_i(x_i-\mu {)}^2}{\nu +n}$
Normal$(\mu ,{\sigma }^2)$	Normal-Scaled Inverse Gamma$(\lambda ,\nu ,\alpha ,\beta )$	$\frac{\nu \lambda +n\bar{x}}{\nu +n},\nu +n,\alpha +\frac{n}{2},\beta +\frac{1}{2}{\sum }_i(x_i-\bar{x}{)}^2+\frac{\gamma (\bar{x}-\lambda {)}^2}{2(n+\gamma )}$

Computation / Markov Chains

A stochastic process $\left\{X_t:t\in T\right\}$ is a collection of random variables.

The process $\left\{X_n:n\in T\right\}$ is a Markov Chain if

$$\mathbf{Pr}\left(X_n=X|X_0,\dots ,X_{n-1}\right)=\mathbf{Pr}\left(X_n=x|X_{n-1}\right)$$

First rewrite the integral to be evaluated $I={\int }_a^{b}h(x)dx$ as follows

$$I={\int }_a^{b}h(x)dx={\int }_a^{b}w(x)f(x)dx$$

where $w(x)=h(x)(b-a)$ and $f(x)=\frac{1}{b-1}$. Since $f$ is the probability density for a uniform r.v. over $(a,b)$, we can write $I={\mathbb{E}}_f\left[w(X)\right]$ where $X\sim \mathrm{Unif}\left(a,b\right)$. If we generate $X_1,\dots ,X_N\sim \mathrm{Unif}\left(a,b\right)$, by LLN,

$$\hat{I}=\frac{1}{N}\sum _{i=1}^{N}w(X_i)\stackrel{p}{\to }\mathbb{E}\left[w(X)\right]=I$$

The goal is the generate sequences ${\mathbf{\theta }}_1,{\mathbf{\theta }}_2,\dots$ from $f(\mathbf{\theta }|\mathbf{y})$. An MCMC scheme will generate samples from the conditional if

$$P({\mathbf{\theta }}_j|{\mathbf{\theta }}_{j-1})f({\mathbf{\theta }}_{j-1}|\mathbf{y})=P({\mathbf{\theta }}_{j-1}|{\mathbf{\theta }}_j)f({\mathbf{\theta }}_j|\mathbf{y})$$

where $P({\mathbf{\theta }}_i|{\mathbf{\theta }}_j)$ is the pdf of ${\mathbf{\theta }}_i$ given ${\mathbf{\theta }}_j$. The LHS is the joint pdf of ${\mathbf{\theta }}_j,{\mathbf{\theta }}_{j-1}$ from the chain, if ${\mathbf{\theta }}_{j-1}$ is from $f(\mathbf{\theta }|\mathbf{y})$. Integrating RHS over $d{\mathbf{\theta }}_{j-1}$ yields $f({\mathbf{\theta }}_j|\mathbf{y})$, so the result states that given ${\mathbf{\theta }}_{j-1}$ is from the correct posterior distribution, the chain generates ${\mathbf{\theta }}_j$ also from the posterior $f(\mathbf{\theta }|\mathbf{y})$.

Basic idea - turn high dimensional problem into several one-dimensional problems. Suppose $(X,Y)$ has joint density $f_{X,Y}(x,y)$. Suppose it is possible to simulate from conditional distributions $f_{X|Y}(x|y)$ and $f_{Y|X}(y|x)$. Let $(X_0,Y_0)$ be starting values. Assuming we have drawn $(X_0,Y_0),\dots ,(X_i,Y_i)$, we generate $(X_{i+1},Y_{i+1})$ as follows

$X_{n+1}\sim f_{X|Y}(x|Y_n)$

$Y_{n+1}\sim f_{Y|X}(y|X_{n+1})$

$\dots$ for multiple parameters

. Let $y_i\stackrel{\mathrm{iid}}{\sim }\mathcal{N}\left(\mu ,{\sigma }^2\right)$. Define precision $\tau =1/{\sigma }^2$

Likelihood: $\mathsf{f}\left(y|\mu ,\tau \right)\sim {\tau }^{n/2}\exp \left(\frac{1}{2}\tau {\sum }_{i=1}^n(y_i-\mu {)}^2\right)$

(Noninformative) Prior: $\pi \mu ,\tau \sim \tau$

Posterior Distribution

$$\pi \left(\mu ,\tau |y\right)\sim {\tau }^{(n/2)+1}\exp \left(-\frac{1}{2}\tau \sum _{i=1}^n(y_i-\mu {)}^2\right)$$

full conditionals:

$\pi (\mu |\tau ,y)=\mathcal{N}\left(\bar{y},(n\tau {)}^{-1}\right)$

$\pi (\tau |\mu ,y)=\Gamma \left(\frac{n}{2},\frac{1}{n}{\sum }_{i=1}^n(y_i-\mu {)}^2\right)$

However, it is typically impossible to write out or sample from full-conditionals.

Let $q(y|x)$ be an arbitrary, friendly distribution we can sample from. The conditional density $q(y|x)$ is called the proposal distribution. MH creates a sequence of observations $X_0,\dots$ as follows Choose $X_0$ arbitrarily. Suppose we have generated $X_0,X_1,\dots ,X_i$. Generate $X_{i+1}$ as follows

• Generate proposal $Y\sim q(y|X_i)$

• Evaluate $r:=r(X_i,Y)$ where

$$r(x,y)=\min \left\{\frac{f(y)}{f(x)}\frac{q(x|y)}{q(y|x)},1\right\}$$

• Set

$$X_{i+1}=\{\begin{array}{ll}Y& \text{w.p}r\\ X_i& \text{w.p}1-r\end{array}$$

Let ${\mathbf{x}}_i$ be observed and ${\mathbf{z}}_i$ missing. The goal is to maximise the log-likelihood of the observed data

$$\ell (\mathbf{\theta })=\sum _{i=1}^n\log p({\mathbf{x}}_i|\theta )=\sum _{i=1}^n\left[\sum _{{\mathbf{z}}_i}p({\mathbf{x}}_i,{\mathbf{z}}_i|\mathbf{\theta })\right]$$

Cannot push log inside the sum because of unobserved variables. EM tackles the problem as follows. Define complete data log likelihood as ${\ell }_c(\mathbf{\theta }):={\sum }_{i=1}^n\log p({\mathbf{x}}_i,{\mathbf{z}}_i|\mathbf{\theta })$. This cannot be computed, since ${\mathbf{z}}_i$ is unknown.

Instead, define $Q(\mathbf{\theta },{\mathbf{\theta }}^{t-1})=\mathbb{E}\left[{\ell }_c(\mathbf{\theta }|\mathcal{D},{\mathbf{\theta }}^{t-1})\right]$ where $t$ is the iteration number and $Q$ is called the auxiliary function.

Expectation (E) Step: Compute $Q(\mathbf{\theta },{\mathbf{\theta }}^{t-1})$, which is an expectation wrt old params ${\mathbf{\theta }}^{t-1}$.

Maximisation (M) Step: Optimise the $Q$ function wrt $\mathbf{\theta }$.

Compute ${\mathbf{\theta }}^t={\mathrm{argmax}}_{\mathbf{\theta }}Q(\mathbf{\theta },{\mathbf{\theta }}^{t-1})$. For MAP estimation, ${\mathbf{\theta }}^t={\mathrm{argmax}}_{\mathbf{\theta }}Q(\mathbf{\theta },{\mathbf{\theta }}^{t-1})+\log p(\mathbf{\theta })$

Probit has the form $p(y_i=1|z_i)=1_{z_i>0}$ where $z_i\sim \mathcal{N}\left({\mathbf{x}}_i^{\prime }\mathbf{\beta },1\right)$ is the latent variable. The complete data log likelihood, assuming a $\mathcal{N}\left(0,{\mathbf{\Sigma }}_0\right)$ prior on $\mathbf{\beta }$.

$$\begin{array}{rl}\ell (\mathbf{z},\mathbf{\beta }|{\mathbf{\Sigma }}_0)& =\log p(\mathbf{y}|\mathbf{z})+\log \mathcal{N}\left(\mathbf{z}|\mathbf{X}\mathbf{\beta },\mathbf{I}\right)+\log \mathcal{N}\left(\mathbf{\beta }|0,{\mathbf{\Sigma }}_0\right)\\ & =\sum _{i=1}^n\log p(y_i|z_i)-\frac{1}{2}(\mathbf{z}-\mathbf{X}\mathbf{\beta }{)}^{\prime }(\mathbf{z}-\mathbf{X}\mathbf{\beta })-\frac{1}{2}{\mathbf{\beta }}^{\prime }{\left({\mathbf{\Sigma }}_0\right)}^{-1}\mathbf{\beta }+\text{const}\end{array}$$

The posterior in the E step is a truncated Gaussian

$$p(z_i|{\mathbf{x}}_i,\mathbf{\beta })=\{\begin{array}{ll}{\mu }_i+\frac{\varphi ({\mu }_i)}{\Phi ({\mu }_i)}& \text{ if}{y}_i=1\\ {\mu }_i-\frac{\varphi ({\mu }_i)}{\Phi ({\mu }_i)}& \text{ if}{y}_i=0\end{array}$$

where ${\mu }_i={\mathbf{x}}_i^{\prime }\mathbf{\beta }$. In the $M$ step, we estimate $\mathbf{\beta }$ using ridge, where $\mathbf{\mu }=\mathbb{E}\left[\mathbf{z}\right]$.

$$\hat{\beta }={\left({\left({\mathbf{\Sigma }}_0\right)}^{-1}+{\mathbf{X}}^{\prime }\mathbf{X}\right)}^{-1}{\mathbf{X}}^{\prime }\mathbf{\mu }$$

Hierarchical Models

Parameters in a prior are modeled as having a distribution that depends on hyperparameters. This results in joint posteriors of the form

$$f(\theta ,\tau |y)\propto \underset{\text{likelihood}}{\underbrace{\mathcal{L}(y|\theta )}}\underset{\text{parameter prior}}{\underbrace{f(\theta |\tau )}}\underset{\text{hyperparameter prior}}{\underbrace{f(\tau )}}$$

Represented by the graphical model $\mathbf{\tau }\to \mathbf{\theta }\to \mathcal{D}$.

We are typically interested in the marginal posterior of $\theta$, which is obtained by integrating the joint posterior w.r.t $\tau$.

By treating $\mathbf{\tau }$ as a latent variable, we allow data-poor observations to borrow strength from data rich ones.

Empirical Bayes

In hierarchical models, we need to compute the posterior on multiple layers of latent variables. For example, for a two-level model, we need

$$p(\mathbf{\eta },\mathbf{\theta }|\mathcal{D})\propto p(\mathcal{D}|\mathbf{\theta })p(\mathbf{\theta }|\mathbf{\eta })p(\mathbf{\eta })$$

We can employ a computational shortcut by approximating the posterior on the hyper-parameters with a point-estimate $p(\mathbf{\eta }|\mathcal{D})\approx {\delta }_{\hat{\mathbf{\eta }}}(\mathbf{\eta })$, where $\hat{\mathbf{\eta }}=\mathrm{argmax}p(\mathbf{\eta }|\mathcal{D})$. Since $\mathbf{\eta }$ is usually much smaller than $\mathbf{\theta }$ in dimensionality, we can safely use a uniform prior on $\mathbf{\eta }$. Then, the estimate becomes

$$\hat{\eta }=\mathrm{argmax}p(\mathcal{D}|\mathbf{\eta })=\mathrm{argmax}\underset{\text{marginal likelihood}}{\underbrace{\left[\int p(\mathcal{D}|\mathbf{\theta })p(\mathbf{\theta }|\mathbf{\eta })d\mathbf{\theta }\right]}}$$

This violates the principle that the prior should be chosen independently of the data, but is a cheap computational trick. This produces a hierarchy of Bayesian methods in increasing order of the number of integrals performed.

Suppose we measure the number of people in various cities $N_i$ and the number of people who died of cancer $x_i$. We assume $x_i\sim \mathrm{Binomial}[N_i,{\theta }_i]$ and want to estimate the cancer rates ${\theta }_i$. The MLE solution would be to either estimate them all separately, or estimate a single $\theta$ for all cities.

The hierarchical approach is to model ${\theta }_i\sim \mathrm{Beta}[a,b]$, and write a joint distribution

$$\begin{array}{rl}p(\mathcal{D},\mathbf{\theta },\mathbf{\eta }|\mathbf{N})& =p(\mathbf{\eta })\prod _{i=1}^n\mathrm{Binomial}[x_i|N_i,{\theta }_i]\mathrm{Beta}[{\theta }_i,\mathbf{\eta }]\\ & \text{analytically integrate out }{\theta }_i\\ & =\prod _{i=1}^n\int \mathrm{Binomial}[x_i|N_i,{\theta }_i]\mathrm{Beta}[{\theta }_i,\mathbf{\eta }]d{\theta }_i\\ & =\prod _{i=1}^n\frac{\mathrm{Beta}[a+x_i,b+N_i-x_i]}{\mathrm{Beta}[a,b]}\end{array}$$

where $\eta :=(a,b)$. We can also put covariates on ${\theta }_i=f({\mathbf{x}}_i^{\prime }\beta )$.

Hierarchy of Bayesianity

Method	Definition
Maximum Likelihood	$\hat{\mathbf{\theta }}={\mathrm{argmax}}_{\mathbf{\theta }}p(\mathcal{D}\mid \mathbf{\theta })$
MAP Estimation	$\hat{\mathbf{\theta }}={\mathrm{argmax}}_{\mathbf{\theta }}p(\mathcal{D}\mid \mathbf{\theta })p(\mathbf{\theta }\mid \mathbf{\eta })$
Empirical Bayes	$\hat{\mathbf{\eta }}={\mathrm{argmax}}_{\mathbf{\eta }}\int p(\mathcal{D}\mid \mathbf{\theta })p(\mathbf{\theta }\mid \mathbf{\eta })d\mathbf{\theta }={\mathrm{argmax}}_{\mathbf{\eta }}p(\mathcal{D}\mid \mathbf{\eta })$
MAP-II	$\hat{\mathbf{\eta }}={\mathrm{argmax}}_{\mathbf{\eta }}\int p(\mathcal{D}\mid \mathbf{\theta })p(\mathbf{\theta }\mid \mathbf{\eta })p(\mathbf{\eta })d\mathbf{\theta }={\mathrm{argmax}}_{\mathbf{\eta }}p(\mathcal{D}\mid \mathbf{\eta })p(\mathbf{\eta })$
Full Bayes	$p(\mathbf{\theta },\mathbf{\eta }\mid \mathcal{D})\propto p(\mathcal{D}\mid \mathbf{\theta })p(\mathbf{\theta }\mid \mathbf{\eta })p(\mathbf{\eta })$

Graphical Models

Any joint distribution can be represented as follows

$$p(x_{1:v})=p(x_1)p(x_2|x_1)p(x_3|x_2,x_1)\dots p(x_v|x_{1:V-1})$$

where $V$ is the number of variables [and we have dropped the parameter vector $\mathbf{\theta }$]. It follows that

The joint distribution $p(\mathbf{x})=p(x_1,\dots ,x_K)$ can be written as

$$p(\mathbf{x})=\prod _{k=1}^{K}p(x_k|{\text{Pa}}_k)$$

where Pa${}_k$ denotes the parent nodes of $x_k$, which are nodes that have arrows pointing to $x_k$.

$X$ and $Y$ are said to be conditionally independent iff the conditional joint can be written as the product of the conditional marginal.

$$X\perp \perp Y|Z⟺p(X,Y|Z)=p(X|Z)p(Y|Z)$$

If the condition $\mathcal{A}\perp \perp \mathcal{B}|\mathcal{C}$, it must be the case that all paths are blocked. All paths are blocked iff

Arrows on the path meet either head to tail or tail to tail at the node, and the note is in the set $\mathcal{C}$

Arrows meet head to head at the node, and neither the node nor any of its descendants is in the set $\mathcal{C}$

Empirical Bayes

We suppose that each player’s MLE value $p_i$ (his batting average in the first 90 tries) is a binomial proportion,

$$p_i\sim \mathrm{Bi}\left(90,P_i\right)/90$$

Here $P_i$ is his true average, how he would perform over an infinite number of tries; TRUTH ${}_i$ is itself a binomial proportion, taken over an average of 370 more tries per player.

At this point there are two ways to proceed. The simplest uses a normal approximation to (7.17)

$$p_i\dot{\sim }\mathcal{N}\left(P_i,{\sigma }_0^2\right)$$

where ${\sigma }_0^2$ is the binomial variance

$${\sigma }_0^2=\bar{p}(1-\bar{p})/90$$

with $\bar{p}=0.254$ the average of the $p_i$ ‘s. Letting $x_i=p_i/{\sigma }_0$, applying (7.13), and transforming back to ${\hat{p}}_i^{\mathrm{JS}}={\sigma }_0{\hat{\mu }}_i^{\mathrm{JS}},$ gives James-Stein estimates

$${\hat{p}}_i^{\mathrm{JS}}=\bar{p}+\left[1-\frac{(N-3){\sigma }_0^2}{\sum {\left(p_i-\bar{p}\right)}^2}\right]\left(p_i-\bar{p}\right)$$

A second approach begins with the arcsin transformation

$$x_i=2(n+0.5{)}^{1/2}{\sin }^{-1}\left[{\left(\frac{n{p}_i+0.375}{n+0.75}\right)}^{1/2}\right]$$