Computing the marginal MAP estimate

Question

Suppose I have some IID Gaussian data with priors for mean and standard deviation:

\begin{align} P(x|\mu,\sigma)&=\prod_{i=1}^n\frac{1}{\sqrt{2\pi\sigma^2}}e^\frac{-(x_i-\mu)^2}{2\sigma^2}\\ p(\mu) &= e^{-\mu},\hspace{3mm}\mu\in(0,\infty)\\ p(\sigma) &= \frac{1}{\sigma^2},\hspace{3mm}\sigma\in(1,\infty) \end{align}

Then suppose I want to evaluate the marginal MAP estimate:

$$\max_{\mu}p(\mu|x)$$

Then I need to compute said marginal distribution

\begin{align} p(\mu|x)\propto p(x|\mu)p(\mu)&=\int p(x|\mu,\sigma)p(\mu)p(\sigma)d\sigma\\ &=\frac{e^{-\mu}}{(2\pi)^\frac{n}{2}}\int_1^\infty\frac{1}{\sigma^{n+2}}\prod_{i=1}^ne^\frac{-(x_i-\mu)^2}{2\sigma^2}d\sigma \end{align}

Supposing that this integral is intractable to compute analytically, I instead apply Laplace's method:

\begin{align} \frac{e^{-\mu}}{(2\pi)^\frac{n}{2}}\int_1^\infty\frac{1}{\sigma^{n+2}}\prod_{i=1}^ne^\frac{-(x_i-\mu)^2}{2\sigma^2}d\sigma&=\frac{e^{-\mu}}{(2\pi)^\frac{n}{2}}\int_1^\infty e^{(n+2)(-\log\sigma-\frac{1}{2(n+2)\sigma^2}\sum_{i=1}^n(x_i-\mu)^2)}d\sigma\\ &=\frac{e^{-\mu}}{(2\pi)^\frac{n}{2}}\int_1^\infty e^{(n+2)f(\sigma)}d\sigma \end{align}

where

$$f(\sigma)=-\log\sigma-\frac{1}{2(n+2)\sigma^2}\sum_{i=1}^n(x_i-\mu)^2$$

The minimum value of $f(\sigma)$ is given by

$$\sigma_0=\sqrt{\frac{1}{n+2}\sum_{i=1}^n(x_i-\mu)^2}$$

and furthermore,

\begin{align} f(\sigma_0) &= -\log\sigma_0 - \frac{1}{2}\\ f''(\sigma_0) &= \frac{-2}{\sigma_0^2} \end{align}

Thus by Laplace's method:

\begin{align} \frac{e^{-\mu}}{(2\pi)^\frac{n}{2}}\int_1^\infty e^{(n+2)f(\sigma)}d\sigma&\approx\frac{e^{-\mu}}{(2\pi)^\frac{n}{2}}\sqrt{\frac{2\pi}{(n+2)\left|f''(\sigma_0)\right|}}e^{(n+2)f(\sigma_0)}\\ &=\frac{e^{-\mu}}{(2\pi)^\frac{n}{2}}\sqrt{\frac{\pi\sigma_0^2}{(n+2)}}e^{-(n+2)(\log\sigma_0+\frac{1}{2})} \end{align}

Optimizing this last expression for $\mu$, we obtain:

$$\hat{\mu}=\frac{2\sum_{i=1}^nx_i-n^2±\sqrt{(n^2-2\sum_{i=1}^nx_i)^2-4n(\sum_{i=1}^nx_i^2-n\sum_{i=1}^nx_i)}}{2n}$$

My point with all this is I would now like to use some other method to verify the correctness of Laplace's method, but I have had trouble finding relevant information on how to use MCMC or some EM, VB or message passing variant for calculating marginal MAP estimates. Can anyone show me another way to calculate this?

Answer 1

I understand the question is asking how to validate via an alternate method that Laplace approximation to the marginalized posterior $p(\mu\vert x)$. The post suggests a couple of methods MCMC, EM, or VB. While these might all be possible the most straightforward method to me appears to be full blown MCMC where you work with the joint distribution $p(\mu, \sigma \vert x)$. To sample from this we write out the full join posterior distribution and then identify the full conditionals for $\mu \vert \sigma, x$ and $\sigma \vert \mu, x$

The join posterior is

$p(\mu, \sigma \vert x ) \propto \frac{\exp[-\mu]}{(2\pi)^{n/2}\sigma^{2+n} }\times \exp[-\frac{-\sum_{i=1}^n(x_i -\mu)^2}{2\sigma^2}] $

now the support of $\mu > 0$ and $\sigma>1$. So the posterior for $\mu \vert \sigma, x$ will be a normal distribution truncated to the positive halfline. The distribution of $\sigma \vert \mu , x$ has support on $\sigma > 1$ and will also be a truncated distribution. This full conditional distribution of $\sigma$ will be a truncated inverse gamma density. All that remains is to determine the parameters of each of the two distributions. We can use the breakdown $\sum_{i=1}^n(x_i -\mu)^2 = (n-1)s^2+n(\mu - \bar{x})^2$ and for the full conditional on $\mu$ we can ignore the first term because this is constant with respect to $\mu$, now to incorporate the $-\mu$ term we expand the square and look only at the terms with $\mu$. The squared term is $\mu^2 n$ and the coefficient on $\mu$ is $\sigma^2 + n\bar{x}$ so this tells us that variance of the full conditional is $\sigma^2/n$ and the full conditional mean divided by the variance of the full conditional is $\sigma^2+n\bar{x}$ so solving for the mean yields $\sigma^4 + \sigma^2\bar{x}$. These are the parameters of the normal distribution truncated to the positive real line.

Similarly the inverse gamma has parameters $\alpha =n+1$ and $\beta = \frac{n}{2}(\bar{x}-\mu)$ here the inverse gamma is truncated to $(1, \infty)$. So the algorithm to sample from the joint posterior is:

For $i$ in $1, \dots, m$

$\sigma_i \sim Inv-Gamma_{[1,\infty)}\left(n+1, \frac{n}{2}(\bar{x}-\mu) \right)$

$\mu_i \sim Normal_{(0, \infty)}\left(\sigma^4 +\sigma^2\bar{x}, \frac{\sigma^2}{n}\right)$

Using the drawn samples from $\mu_i$ and a kernal density estimator, you should be able to determine if you are in the correct vicinity of the posterior with your Laplace approximation. Now for simulating random variates from a truncated distribution there are a couple of standard methods, here is one approach that will work numerically. The R programming language has standard cumulative densities for Inverse Gammas and Normal distributions as well as the quantile/inverse cdf functions.