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

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.

• This is an interesting setup; I'm curious as to how this prior, nonsingular as it is, came about, and I wonder if you could alter it to make it full rank with the "extended" part of $J$ formed to be a relatively uninformative part of the covariance matrix. Commented Aug 16, 2023 at 15:20
• Since you use a normal distribution for the parameter $\theta$, in what way or how is it improper? Commented Aug 16, 2023 at 18:00
• Or you mean that you use a prior with is $k$ time a gaussian kernel, with $k$ such that the integral of the product of $k$ with the Gaussian kernel is infinite? Commented Aug 16, 2023 at 18:06
• @Fiodor1234 - since $A$ is singular, the prior is actually on a collection of linear combinations of the $\theta_i$ that is not $1-1$ with respect to $\theta$; for example, imagine that the prior is on the sum of the $\theta_i$, not the $\theta_i$ themselves. Commented Aug 16, 2023 at 19:07
• @SextusEmpiricus - see my example in the comment above or the explanation by jbowman. I believe there are endless examples where the prior knowledge only applies to a subset of the model we want to estimate. If the model is a FIR and the prior information is that of the model gain, we have that the prior is one on the sum of the parameters, as jbowman exemplifies. It is improper since the prior covariance matrix on $\theta$ is singular. Commented Aug 17, 2023 at 11:15

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.

• No, this is not where my confusion stems from, and not why division by zero occurs. Computing the posterior $p(\mu | y_1)$ is perfectly possible as with "any" problem in Bayesian methods, and in my case this posterior is proper by assumption. Why do you claim that $f(y_1)$ is improper? Commented Aug 18, 2023 at 7:47
• @smallStackBigFlow because $f(\mu)$ is improper and due to symmetry we can say that $f(y_i) = f(\mu)$. Commented Aug 18, 2023 at 9:17
• Or work out the integral $$f(y_i) = \int f(y_i | \mu) f(\mu) d\mu = f(\mu) \underbrace{ \int f(y_i | \mu) d\mu }_{=1} = f(\mu)$$ the term $f(\mu)$ can be moved outside the integral because it is constant. Commented Aug 18, 2023 at 9:24
• Also note that I didn't state that the posterior $p(\mu|y_1)$ or $f(\mu|y_1)$ is improper, which your comment suggests. It is $f(y_1)$ that is improper. Commented Aug 18, 2023 at 9:53
• @smallStackBigFlow I have updated with a slightly less simple problem. What I am considering is that the problem with the matrix $A$ can be seen as effectively a part that is improper and a part that is proper. E.g., if we have a prior for only the sum $\theta_1+ \theta_2$ and then we could consider a reparameterization with $\tau_1 = \theta_1 + \theta_2$ and $\tau_2 = \theta_1 - \theta_2$, and the problem is like a separate treatment of $\tau_1$ with a proper prior, and $\tau_2$ with an improper prior. As an example for the problem it might be easier to start right away with these variables. Commented Aug 18, 2023 at 13:25