Note: Marginal likelihood for Bayesian models with a Gaussian-approximated posterior

I first learned this solution from Botond Cseke . I'm not sure where it originates; It is essentially Laplace's method for approximating integrals using a Gaussian distribution, where the parameters of the Gaussian distribution might come from any number of various approximate inference approaches.

If I have a Bayesian statistical model with hyperparameters $\mathrm{\Theta }$, with a no closed-form posterior, how can I optimize $\mathrm{\Theta }$?

Consider a Bayesian statistical model with observed data $\mathbf{y}\in \mathcal{Y}$ and hidden (latent) variables $\mathbf{z}\in \mathcal{Z}$, which we infer. We have a prior on $\mathbf{z}$, $Pr(\mathbf{z};\mathrm{\Theta })$, and a model for the probability of $\mathbf{y}$ given $\mathbf{z}$ (likelihood), $Pr(\mathbf{y}|\mathbf{z};\mathrm{\Theta })$. The prior and likelihood are controlled by "hyperparameters" $\mathrm{\Theta }$, which we would like to estimate. Recall that Bayes theorem states:

$$\begin{array}{}\text{(1)}& \begin{array}{r}Pr(\mathbf{z}|\mathbf{y};\mathrm{\Theta })=Pr(\mathbf{y}|\mathbf{z};\mathrm{\Theta })\frac{Pr(\mathbf{z};\mathrm{\Theta })}{Pr(\mathbf{y};\mathrm{\Theta })}\end{array}\end{array}$$

It is common for the posterior $Pr(\mathbf{z}|\mathbf{y};\mathrm{\Theta })$ to lack a closed-form solution. In this case, one typically approximates the posterior with a more tractable distribution $Q(\mathbf{z})\approx Pr(\mathbf{z}|\mathbf{y};\mathrm{\Theta })$. Common ways of estimating $(\mathbf{z})$ include the Laplace approximation , variational Bayes , expectation propagation , and expectation maximization algorithms. The only approximating distribution in common use for high-dimensional $\mathbf{z}$ is the multivariate Gaussian (or some nonlinear transformation thereof), which succinctly captures joint statistics with limited computational overhead. Assume we have an inference procedure which returns the approximate posterior $Q(\mathbf{z})=\mathcal{N}({\mathit{\mu }}_q,{\mathbf{\Sigma }}_q)$.

We optimize the hyperparameters "$\mathrm{\Theta }$" of the prior kernel to maximize the marginal likelihood of the observations $\mathbf{y}$

$$\begin{array}{}\text{(2)}& \begin{array}{rl}\theta & \leftarrow \underset{\mathrm{\Theta }}{\mathrm{argmax}}Pr(\mathbf{y};\mathrm{\Theta })\\ Pr(\mathbf{y};\mathrm{\Theta })& ={\int }_{\mathcal{Z}}Pr(\mathbf{y},\mathbf{z};\mathrm{\Theta })d\mathbf{z}={\int }_{\mathcal{Z}}Pr(\mathbf{y}|\mathbf{z})Pr(\mathbf{z};\mathrm{\Theta })d\mathbf{z}\end{array}\end{array}$$

Except in rare special cases, this integral does not have a closed form. However, we have already obtained a Gaussian approximation to the posterior distribution, $Q(\mathbf{z})\approx Pr(\mathbf{z}|\mathbf{y};\mathrm{\Theta })$. If we replace $Pr(\mathbf{z}|\mathbf{y};\mathrm{\Theta })$ with our approximation $Q(\mathbf{z})$ in this equation, we can solve for (an approximation) of $Pr(\mathbf{y};\mathrm{\Theta })$:

$$\begin{array}{}\text{(3)}& \begin{array}{rl}Q(\mathbf{z})\approx Pr(\mathbf{y}|\mathbf{z})\frac{Pr(\mathbf{z};\mathrm{\Theta })}{Pr(\mathbf{y};\mathrm{\Theta })}& \Rightarrow Pr(\mathbf{y};\mathrm{\Theta })\approx Pr(\mathbf{y}|\mathbf{z})\frac{Pr(\mathbf{z};\mathrm{\Theta })}{Q(\mathbf{z})}\end{array}\end{array}$$

Working in log-probability, and evaluating the expression at the (approximated) posterior mean $\mathbf{z}={\mathit{\mu }}_q$, we get

$$\begin{array}{}\text{(4)}& \begin{array}{rl}\ln Pr(\mathbf{z}={\mathit{\mu }}_q;\mathrm{\Theta })& =-\frac{1}{2}\{\ln |2\pi {\mathbf{\Sigma }}_z|+({\mathit{\mu }}_q-{\mathit{\mu }}_z{)}^{\mathrm{\top }}{\mathbf{\Sigma }}_z^{-1}({\mathit{\mu }}_q-{\mathit{\mu }}_z)\}\\ \ln Q(\mathbf{z}={\mathit{\mu }}_q)& =-\frac{1}{2}\{\ln |2\pi {\mathbf{\Sigma }}_q|+({\mathit{\mu }}_q-{\mathit{\mu }}_q{)}^{\mathrm{\top }}{\mathbf{\Sigma }}_q^{-1}({\mathit{\mu }}_q-{\mathit{\mu }}_q)\}=-\frac{1}{2}\ln |2\pi {\mathbf{\Sigma }}_q|\\ \\ \ln Pr(\mathbf{y};\mathrm{\Theta })& \approx \ln Pr(\mathbf{y}|\mathbf{z}={\mathit{\mu }}_q;\mathrm{\Theta })+\ln Pr(\mathbf{z}={\mathit{\mu }}_q;\mathrm{\Theta })-\ln Q(\mathbf{z}={\mathit{\mu }}_q)\\ & =\ln Pr(\mathbf{y}|\mathbf{z}={\mathit{\mu }}_q;\mathrm{\Theta })-\frac{1}{2}\{\ln |{\mathbf{\Sigma }}_q^{-1}{\mathbf{\Sigma }}_z|+({\mathit{\mu }}_q-{\mathit{\mu }}_z{)}^{\mathrm{\top }}{\mathbf{\Sigma }}_z^{-1}({\mathit{\mu }}_q-{\mathit{\mu }}_z)\}.\end{array}\end{array}$$

This is quite tractable to compute.