Hierarchical gaussian models, empirical bayes and shrinkage

This article is the second of the Bayesian series. Previous episode here.

You probably already know the Bayes theorem or stumble upon Naive Bayes when comparing machine learning models. Maybe, like me, you feel that you barely saw the shores of the continent that is the Bayesian land and want to set foot in it.

I wrote this article from a Stanford course with the intent for you to understand the relationships between all Bayesian concepts without having to go through an entire classroom (like I did). As a matter of readability, I hide most of the maths in dropdowns. This way, you can follow the logic easily and come back to deep dive into formulas later.

So far we have dealt with Bayesian analysis of very simple models, mainly using conjugate priors, with few parameters to estimate. In the coin tossing example, we had just one parameter, the probability of heads. In the Gaussian model, we had two parameters, the mean $μ$ and variance $σ^2$. However, in modern applications one typically must estimate many parameters, using more complicated hierarchical models.

Suppose that we have $N$ observations each from independent Gaussians all with same variance 1 but with different means

$y_i \;\sim^{ind}N(\mu_i,1)$​

This is a high-dimensional model: there are as many parameters as there are data. Another way to write it with N-vectors:

$y\;\sim\;N(I\mu,I)$​

$R_{frequentist}-R_{Bayesian}=\frac{N}{\tau^2+1}$, which can be very large if τ is small.

Obviously we want estimators that have low risk. For any loss function $L$ we know that:

Where ξ are prior hyperparameters. However frequentists don’t use prior, so instead we can use admissibility and minimax.

An estimator $\delta^*$ is inadmissible if there exists another estimator $\delta$ for which:

$R(\theta,\delta)\leq R(\theta,\delta^*)$​

So an estimator is admissible if it is not inadmissible. $\delta$ dominate $\delta^*$ if $R(θ, δ)≤R(θ, δ^*)$ for every $\theta$.

If we don’t know how to choose the parameter $\tau$ of $\mu$, we try to estimate $\mu$ using empirical Bayes via the James-Stein estimator.

$\hat{\mu}^{(JS)}=(1-\frac{N-2}{||y||^2})\;y$​

So why one would use the MLE is high-dimension? Because:

In high dimension, suppose we generate values from

$y_i\;\sim\;N(\mu_i,1)\;for\;i=1,...,11$​

We put

$\mu_i=(5,-1,-0.75,...,0.75,1)$​

For each simulation replicate, we estimate μ by the MLE and by JS. We find:

We could conclude at this point that one should always use empirical Bayes to choose a prior. However, when it comes to estimating a parameter, the Bayesian approach is to put a prior on it.

Thus if our problem is

$y\;\sim\;N(\mu,I)\\\mu\;\sim\;N(0,\tau^2I)\\\tau^2\;\sim\;p(\tau^2)$​

now our problem is to estimate $p(\tau^2)$. Since we want to estimate $\tau^2$ from our data, we have to place an objective prior on $\tau^2$: $p(τ^2)=1$​

We need some alternative way to do integration that doesn’t involve actually computing all of the integrals analytically.

We still have the same problem. Now we know that getting the posterior is hard if we can’t find the normalising constant through some known distribution.

Another choice is the half Cauchy prior on τ

$p(\tau)\propto\frac{1}{1+\tau^2}$​

We choose to have an Inverse Gamma prior on $\tau^2$ and a Normal prior on $m$.

However, in high-dimensions we don’t really care about the whole posterior, but more about some functional of it, like the posterior expectations of a function f:

$\int\frac{f(\theta)p(y|\theta)p(\theta)}{p(y)}d\theta$​

It turns out this integral can be approximated using Markov chain Monte Carlo (MCMC). It consists in sampling the desired distribution by recording states from the Markov chain, which has the desired distribution as its equilibrium distribution.

A Markov chain is a sequence of random variable ignoring the whole history (independent) except for the last event:

$p(X_N|X_0,...,X_{N-1})=p(X_N|X_{N-1})$​

We hop from one state $n$ to state $n+1$ via a transition kernel $K$. Some $K$ have invariant distribution.

We want to construct $K$ for arbitrary distributions. We need to compute the posterior without knowing the normalising constant $p(y)$.

In Gibbs sampling the idea is to break the problem of sampling from the high-dimensional joint distribution into a series of samples from low-dimensional conditional distributions.

A nice implementation can be found here:

If one of the Gibbs conditionals is not conjugate, a Metropolis-Hasting step is used in the Gibbs iteration.

Here is an implementation: MCMC sampling for dummies

This new step ensure the transition kernel is reversible, which guarantee that the Markov chain converges to the correct invariant measure.

The core idea is to compute an acceptance ratio, which is the probability to jump to the next step.

We can diagnosis convergence via:

We plot the autocorrelation of the chain, ideally quickly decreasing to zero (the chain itself looks like white noise).

However the shortcomings of this idea is to only consider the effect of auto-correlations for the identity function.

We need to bound autocorrelations for a huge class of functions if we want to converge in total variation

We see how Bayesian linear regression

$z_{N×1}=W_{N×d}\beta_{d×1}+\epsilon_{N×1}\\with\;\epsilon\sim N(0,\sigma^2I)$​

easily leads to ridge estimate

$\hat{\beta}_{ridge}=\mathbb{E}[\beta|W,\sigma^2,z,\tau^2]=\Big(W^TW+\frac{1}{\tau^2}I\Big)^{-1}W^Tz$​

where z is the response, W is the matrix, β is the vector of regression coefficients, and ϵ is the vector of random error.

More generally we have Tikhonov regularization

$\hat{\beta}_{ridge}=(W^TW+A^{-1})^{-1}W^Tz$​

where A symmetric, positive definite.

We regularise because:

Another option to solve Bayesian linear regression is using a conjugate prior.