Laplace approximation

Section of Model selection for stochastic models, part of Statistical learning

Type	Section
Taught in	Model selection for stochastic models
Part	Statistical learning
Questions using it	1

The integral is what has to be approximated. Expand \log\mathbb{P}(X_{\mathrm{all}},\Psi_k\mid M_k) to second order about a mode \widetilde\Psi_k. The gradient vanishes at a mode, so the linear term is absent and the quadratic term is the negative Hessian:

\log\mathbb{P}(X_{\mathrm{all}},\Psi_k\mid M_k) \approx \log\mathbb{P}(X_{\mathrm{all}},\widetilde\Psi_k\mid M_k) -\frac12(\Psi_k-\widetilde\Psi_k)^\top H(\widetilde\Psi_k)\,(\Psi_k-\widetilde\Psi_k).

Here H is the negative Hessian of the log joint density. An interior mode guarantees a positive semidefinite negative Hessian. Strict positive definiteness is an additional assumption. The function -x^4 has a maximum at zero and zero second derivative there. Laplace’s Gaussian approximation needs nonzero curvature in every parameter direction. A singular direction breaks this particular argument.

Substitute the quadratic approximation into the integral for k^*:

k^*\approx\arg\max_k \mathbb{P}(M_k)\,\mathbb{P}(X_{\mathrm{all}},\widetilde\Psi_k\mid M_k) \int_{\Theta_k} \exp\Bigl\{-\tfrac12(\Psi_k-\widetilde\Psi_k)^\top H(\widetilde\Psi_k)\,(\Psi_k-\widetilde\Psi_k)\Bigr\} \,d\Psi_k.

Up to the normalising constant, the integrand is the density of the d_k-dimensional normal distribution with mean \widetilde\Psi_k and covariance H(\widetilde\Psi_k)^{-1}. The domain \Theta_k is not all of \mathbb{R}^{d_k}. Weights lie in (0,1), variances are positive, and covariance matrices are constrained, so the integral of that normal kernel over the legal region is not exactly 1.

The normalising factor of the kernel is (2\pi)^{d_k/2}|H|^{-1/2}. Extending the local integral to all of \mathbb{R}^{d_k} is justified asymptotically when the mass concentrates near an interior, identifiable parameter and the information per observation approaches a positive definite matrix. For a smooth positive prior, the posterior mode and the maximum likelihood estimator are then asymptotically close. They are not exactly equal. These assumptions keep the boundary away from the shrinking region that contributes most of the integral (Drton and Plummer, 2017). Use |A^{-1}|=|A|^{-1} to write the approximate integral as

\int_{\Theta_k} \exp\Bigl\{-\tfrac12(\Psi_k-\widetilde\Psi_k)^\top H(\widetilde\Psi_k)\,(\Psi_k-\widetilde\Psi_k)\Bigr\}\,d\Psi_k \approx (2\pi)^{d_k/2}\,\bigl|H(\widehat\Psi_k)\bigr|^{-1/2}.

The approximate MAP index is therefore

k^*\approx\arg\max_k \mathbb{P}(M_k)\,\mathbb{P}(X_{\mathrm{all}},\widehat\Psi_k\mid M_k)\, (2\pi)^{d_k/2}\,\bigl|H(\widehat\Psi_k)\bigr|^{-1/2}.

The logarithm is strictly increasing, so it preserves the maximiser. Factor the joint probability at \widehat\Psi_k into the likelihood and the prior density:

\begin{align*} k^* &\approx\arg\max_k \log\mathbb{P}(M_k)+\log\mathbb{P}(X_{\mathrm{all}}\mid\widehat\Psi_k,M_k)+\log\mathbb{P}(\widehat\Psi_k\mid M_k)\\ &\qquad+\frac{d_k}{2}\log(2\pi)-\frac12\log\bigl|H(\widehat\Psi_k)\bigr|\\ &\approx\arg\max_k \log\mathbb{P}(M_k)+\ell(\widehat\Psi_k;X_{\mathrm{all}},M_k)+\log\mathbb{P}(\widehat\Psi_k\mid M_k)\\ &\qquad+\frac{d_k}{2}\log(2\pi)-\frac12\log I(\widehat\Psi_k,X_{\mathrm{all}}). \end{align*}

The observed information I is the negative Hessian of the log-likelihood alone. The negative Hessian of the log joint density also contains prior curvature. In this regular large-sample calculation the likelihood contribution is of order n, while a fixed smooth prior contributes order one. Replacing the joint curvature by observed information is therefore another approximation, which is what the second line records. An EM output must also be checked before treating it as the required maximiser.

Used in

Taught in Model selection for stochastic models
Questions that use it (1) — The Laplace integral in one dimension (answer)
All 7 questions on Model selection for stochastic models

Used in

See also