Doing empirical Bayes with improper prior - marginals that do not exist?

Question

I am considering a Bayesian linear model for which the prior is not proper. The model is as usual $y = X \theta + w$ where $w \sim N(0, \sigma^2)$, and $\theta, \sigma^2$ are unknown. The distribution is a normal-inverse Gamma for parameters $\theta, \sigma^2$ (the normal part is improper). There are hyperparameters in the prior that would be nice to tune with empirical Bayes, though I am running into a problem.

I am not able to calculate what $p(y)$ is, i.e. the integral $$p(y) = \int \int p(y | \theta, \sigma^2) p(\theta | \sigma^2) p(\sigma^2) d\theta d\sigma^2$$ divides by zero, and I do not know if this is the right way to do it or where the mistake is. I have followed derivations similar to these to find $p(y)$. If we arrive at the equation we can set the hyperparameters so that we maximize the probability (as described on EB wikipedia).

The prior is constructed from some information that $$ p(\theta | \sigma^2 ) \propto (\sigma^2)^{-q/2} \exp ( - (v - A \theta)^T J (v- A\theta) /(2\sigma^2)) $$ where $v, A, J$ are known, $J$ is positive definite, but $A$ is singular. Therefore this does not define a proper distribution over $\theta$, there "is no" covariance matrix, though there is a singular precision matrix. This is not a problem in order to compute a posterior, we must just assume that $ X^TX + A^T JA $ is invertible, which we assume, then the posterior is proper. The hyperparameters to be tuned are embedded in $J$.

Calculating $p(y)$ as they do in the source leads to dividing by the determinant of the prior covariance matrix, but this is 0, since the matrix is singular. I still think that as long as the posterior is proper, which we assume, we can do empirical Bayes like described, but how? what am I missing? Any resources that consider this scenario?

In the integral, $ p(\theta | \sigma^2) $ is not a proper distribution, since it can not be normalized to 1, so there is an abuse of notation there. This is perhaps partly where the problem comes from, but I still find it hard to believe that empirical Bayes can not be applied when the posterior is proper.

Answer 1

I am not able to calculate what $p(y)$

If the prior of the parameters is improper, then the prior probability of the observations can be improper as well.

So trying to compute the prior probability of the observation $f(y)$ can be problematic and might lead to a division by zero (or infinity).

A simpler model might make this more intuitive.

Consider an observation $Y$ that follows a normal distribution

$$Y \sim N(\mu,\sigma)$$

with fixed $\sigma = 1$ and improper prior for $\mu$

$$f(\mu) \propto 1$$

then given an observation $y$ we have the posterior

$$f(\mu|y) = \frac{1}{\sqrt{2\pi}}e^{-\frac{1}{2}(y-\mu)^2}$$

• One may try to derive this via $$f(\mu|y) = f(y|\mu) \frac{f(\mu)}{f(y)}$$ but here, the ratio $\frac{f(\mu)}{f(y)}$ is difficult to compute because both sides are improper distributions. (And this is probably what you have with the division by zero)

  The distribution $f(y)$ being improper just like $f(\mu)$ follows from: $$f(y) = \int f(y | \mu) f(\mu) d\mu = f(\mu) \underbrace{ \int f(y | \mu) d\mu }_{=1} = f(\mu)$$

• Alternatively we can use $$f(\mu|y) = \frac{f(y|\mu) f(\mu)}{\int_{\forall \mu}f(y|\mu) f(\mu) \,\text{d}\mu}$$

Using the second formulation may be an easier approach and could work as well for your case. Don't compute $f(y)$ and perform the normalization at the end after multiplying all the terms in $f(y|\mu) f(\mu)$.

---

Calculating $p(y)$ as they do in the source leads to dividing by the determinant of the prior covariance matrix, but this is $0$

The problem is not a division by zero, but a division by infinity.

The constant prior $f(\mu) \propto \text{constant}$ can be seen as the limit at infinity, $\sigma_0 \to \infty$, of the conjugate prior for normal distributed observations with known variance

$$f(\mu|\sigma_0) = (2\pi\sigma_0^2)^{-1} \exp(\frac{1}{2}(\mu-\mu_0)^2/\sigma_0^2)$$

In the example from the question this gets complicated further than the simple example above, because the parameter $\mu$ is multivariate and the improper part is in an extra variable.

We could analyze this possibly instead with the following model:

$$Y_1 \sim N(\mu_1,\sigma)\\ Y_2 \sim N(\mu_2,\sigma)$$

with prior

$$f(\mu_1,\mu_2|\sigma_{0,1},\sigma_{0,2}) = \frac{1}{\sqrt{2\pi\sigma_{0,1}}}e^{-\frac{1}{2}(\mu_1)^2/\sigma_{0,1}^2} \frac{1}{\sqrt{2\pi\sigma_{0,2}}}e^{-\frac{1}{2}(\mu_2)^2/\sigma_{0,2}^2}$$

And when, say $\mu_2$ is an improper prior $\sigma_{0,2} \equiv \infty$ and $\sigma_{0,1} =1$ then it is natural to use:

$$f(\mu_1,\mu_2) = \frac{1}{\sqrt{2\pi}}e^{-\frac{1}{2}(\mu_1)^2}$$

Instead, we might also have a treatment of $\mu_i$ together with a single prior parameter $\sigma_0$ (that is fixed for a simple example, but could also follow the inverse gamma distribution)

$$f(\boldsymbol{\mu}|\sigma_{0}) = \frac{1}{\sqrt{2\pi \sigma_0 det(\Sigma)}}e^{-\frac{1}{2}(\boldsymbol{\mu} \Sigma^{-1} \boldsymbol{\mu})}$$

Potentially we could consider $\Sigma^{-1} = \begin{bmatrix} 1 & 0\\0&0\end{bmatrix}$ which relates to the prior on $\mu_2$ being an improper prior.

$$f(\boldsymbol{\mu}|\sigma_{0}) = \frac{1}{\sqrt{2\pi \sigma_0 det(\Sigma)}}e^{-\frac{1}{2}\mu_1^2}$$

I believe that here we have a example situation with a more similar same problem as in the question as there is this division by the $det(\Sigma)$. I guess that we could ignore this division in the formulas as it is effectively just normalisation constant and it becomes a division by zero (or is it actually a division by infinity?) because the $det(\Sigma)$ includes an improper part.