Parameter distribution of $\theta$ from a rectangular matrix multiplication $C\theta$

Question

I am struggeling to see where this problem fits - i.e. what topics this problem relates to, so I am not able to find the right literature. I want to use some particular information as a prior to a model.

In general I want to say something about the distribution of $\theta$ given that I know that $C\cdot \theta \sim \mathcal{N}(\mu, \Sigma)$, when $C$ is rectangular (e.g. it can be a row vector of 1's). I want to use this information as a prior on the model parameters $\theta$. If it helps we can assume that the model is linear in the parameters $Y = X\theta + \epsilon$, with Gaussian noise $\epsilon$, and we have sufficiently many observations of $Y, X$ (I say this since the matrix $C$ may not regularize the problem sufficiently to guarantee a unique solution to the parameter estimate for any number of observations), so we can obtain an estimate of $\theta$ which is Gaussian.

The simplest form of my problem is when $\theta \in \mathbb{R}^2$ and $(1, 1) \cdot \theta \sim \mathcal{N} (\mu, \sigma^2)$ for some given $\mu, \sigma^2$. Then I guess the (prior) distribution on $\theta$ should be a Gaussian along the line $\theta_2 = \mu - \theta_1$. Can we state that somehow? Such a prior should "pull" the parameter estimate towards this subspace defined by $(1,1)$, and the degree of "pulling" is determined by the relative uncertainty of the prior vs the observations of course.

How can I use a prior like the above to compute a posterior? Are there conditions on $C$ (such as full row rank) for this to work? Could you show how the calculations of the posterior are carried out?

Answer 1

In a Bayesian linear regression we can indeed encode the wanted form of relation by considering a prior covariance which is "infinite". This is sometimes called a diffuse prior or a partially diffuse prior.

Consider as in OP the linear regression $\mathbf{y} = \mathbf{X}\boldsymbol{\theta} + \boldsymbol{\varepsilon}$ with a known noise variance $\sigma_{\varepsilon}^2$ where $\mathbf{X}$ is $n \times p$. Remind that the multivariate normal is a conjugate distribution for the parameter $\boldsymbol{\theta}$. If the prior for $\boldsymbol{\theta}$ is the multivariate normal with mean $\mathbf{b}_0$ and covariance $\boldsymbol{\Gamma}_0$ corresponding to a precision matrix $\mathbf{B}_0 = \boldsymbol{\Gamma}_0^{-1}$ then the posterior is also normal with mean $\mathbf{b}_n$ and precision $\mathbf{B}_n$ given by

\begin{align*} \mathbf{B}_n & = \mathbf{B}_ 0 + \sigma^{-2}_\varepsilon \mathbf{X}^\top \mathbf{X}, \\ \mathbf{B}_n \mathbf{b}_n &= \mathbf{B}_0 \mathbf{b}_0 + \sigma^{-2}_\varepsilon \mathbf{X}^\top \mathbf{y}. \end{align*}

The so-called diffuse prior corresponds $\boldsymbol{\Gamma}_0 = \lambda \mathbf{I}_p$ with $\lambda \to \infty$ or $\mathbf{B}_0 = \mathbf{0}$. If $\mathbf{X}$ has full column rank the posterior distribution tends to the normal with covariance $\sigma^2_\varepsilon [\mathbf{X}^\top\mathbf{X}]^{-1}$. This can be regarded as a proper posterior. If $\mathbf{X}$ has rank $< p$, we get an improper posterior which still can be used with some limitations. For instance if all the rows of $\mathbf{X}$ are identical to one vector $\mathbf{x}$ then we have a proper posterior for the corresponding response $\mathbf{x}^\top \boldsymbol{\theta}$ and we can make a proper prediction for a new observation corresponding to this design vector.

A partially diffuse prior can be obtained by using a prior covariance $$ \boldsymbol{\Gamma}_0 = \boldsymbol{\Gamma}_0^{[0]} + \lambda \, \boldsymbol{\Gamma}_0^{[1]}, \qquad \lambda \to \infty \tag{1} $$ where $\boldsymbol{\Gamma}_0^{[0]}$ and $\boldsymbol{\Gamma}_0^{[1]}$ are positive semi-definite matrices. Depending on the design matrix $\mathbf{X}$ and the matrices $\boldsymbol{\Gamma}_0^{[0]}$ and $\boldsymbol{\Gamma}_0^{[1]}$ we can get either an improper posterior or a proper posterior. Even an improper posterior can be used with some limitations. If $\boldsymbol{\Gamma}_0^{[1]}$ has full rank the prior is equivalent to the diffuse prior with covariance $\lambda \mathbf{I}_p$, and the interest of the partially diffuse prior is when $\boldsymbol{\Gamma}_0^{[1]}$ is rank deficient, see example below.

Back to the OP we can consider a prior information $\mathbf{C} \boldsymbol{\theta} \sim \mathcal{N}(\boldsymbol{\mu},\, \boldsymbol{\Sigma})$ for some $r \times p$ matrix $\mathbf{C}$ with full row rank. A special case is in this question corresponding to a zero covariance in the constraint. We can get a (linear) re-parameterization of the linear regression as $$ \mathbf{y} = \mathbf{Z}_1 \boldsymbol{\gamma}_1 + \mathbf{Z}_2 \boldsymbol{\gamma}_2 + \boldsymbol{\varepsilon} $$

where the parameter $\boldsymbol{\gamma}$ stacks $\boldsymbol{\gamma}_1 := \mathbf{C} \boldsymbol{\theta}$ and $\boldsymbol{\gamma}_2$ which is a vector with length $p-r$. We can choose $\boldsymbol{\gamma}_1$ and $\boldsymbol{\gamma}_2$ to be a priori independent with the wanted (proper) prior for $\boldsymbol{\gamma}_1$, and with $\boldsymbol{\gamma}_2$ having the prior mean zero and the prior covariance $\lambda \, \mathbf{I}_{p-r}$ with $\lambda \to \infty$. We can derive the corresponding posterior for $\boldsymbol{\gamma}$ hence for $\boldsymbol{\theta}$. This posterior can be either proper or improper depending on the matrices $\mathbf{X}$ and $\mathbf{C}$. It is easy to see that if the matrix obtained by stacking the rows of $\mathbf{X}$ and those of $\mathbf{C}$ has full rank then the posterior will be proper.

Example Suppose that $p=2$ and that $\mathbf{C}$ is the row matrix $[1,\, 1]$ as in OP so that $\mathbf{C}\boldsymbol{\theta}$ is a scalar r.v. To fix the ideas, we consider the simple linerar regression $y = \theta_1 + \theta_2 x + \varepsilon$ so that $\theta_1 + \theta_2$ is the value for $x= 1$ and $\theta_2 - \theta_1$ is the opposite of the value for $x = -1$. We take $\mu= -2$: if the prior variance $\Sigma$ is small the fitted line should nearly pass by the point $[1, \, -2]$ shown in blue, along with a $\pm 2 \sqrt{\Sigma}$ "error bar". Whatever be $\Sigma$ we should not be informative on the value for $x = -1$. The reparameterization can be $\boldsymbol{\gamma} = \mathbf{G}\boldsymbol{\theta}$ where

$$ \mathbf{G} = \begin{bmatrix} 1 & 1 \\ - 1 & 1 \end{bmatrix}, \qquad \mathbf{G}^{-1} = \frac{1}{2} \, \begin{bmatrix} 1 & -1 \\ 1 & 1 \end{bmatrix}, $$

Let us choose the prior covariance for $\boldsymbol{\gamma}$ as above: $\text{Cov}(\gamma_1) = \Sigma$ and $\text{Cov}(\gamma_2) = \lambda$. The corresponding prior covariance for $\boldsymbol{\theta}$ takes the form (1) above

$$ \text{Cov}(\boldsymbol{\theta}) = \mathbf{G}^{-1} \begin{bmatrix} \Sigma & 0 \\ 0 & \lambda \end{bmatrix} \mathbf{G}^{-\top} = \frac{\Sigma}{4} \, \begin{bmatrix} 1 & 1 \\ 1 & 1 \end{bmatrix} + \frac{\lambda}{4} \begin{bmatrix} 1 & -1 \\ -1 & 1 \end{bmatrix} $$

and note that the matrix coefficient of $\lambda$ has rank one. We can either use a software that accepts a partially diffuse prior, or simply take a large value for $\lambda$ as illustrated here. When $\Sigma$ is large the effect of the prior is negligible, and when $\Sigma$ gets smaller the regression line is driven towards the point $[1, \mu]$.

pdReg <- function(X, y, varNoise = 1,
                  b0 = rep(0, ncol(X)),
                  Gamma0, Gamma1, lambda = 1e4) {
    B0 <- solve(Gamma0 + lambda * Gamma1)
    Bn <- B0 + crossprod(X) / varNoise
    bn <- solve(Bn, B0 %*% b0 + crossprod(X, y) / varNoise)
    list(bn = bn, Bn = Bn)
}
set.seed(123)
## Observations with a specific theta
n <- 10
varNoise <- 0.1
theta <- c(-1, 4)
x <-  seq(-1, 1, length = n)
X <- cbind(Cst = 1, x = x)
y <- X %*% theta + rnorm(n, sd = sqrt(varNoise))
## standard OLS fit
fit0 <- lm(y ~ X - 1)
plot(y ~ x, pch = 16, col = "orangered",
     main = "OLS fit (black) and Bayesian fits (colors)")
abline(coef(fit0))

## Expectation 'mu' and variance (not sd) of C %*% theta
mu <- -2
Sigmas <- c(1e0, 1e-1, 4e-4)
cols <- c("SpringGreen3", "Orchid", "SteelBlue3") 
points(x = 1, y = mu, pch = 16, col = "SteelBlue", cex = 1.5)

for (i in seq_along(Sigmas)) {
    Sigma <- Sigmas[i]
    ## Matrices in the partially diffuse prior
    Gamma0 <- Sigma * matrix(c(1, 1, 1, 1), nrow = 2) / 4
    Gamma1 <-  matrix(c(1, -1, -1, 1), nrow = 2) / 4
    fit <- pdReg(X = X, y = y, varNoise = varNoise,
                 b0 = c(mu, 0), Gamma0 = Gamma0, Gamma = Gamma1, lambda = 1e8)
    abline(fit$bn, col = cols[i], lwd = 2, lty = i)
    segments(x0 = 1 + (i-2) * 0.02, y0 = mu - 2 *sqrt(Sigma), y1 = mu + 2 *sqrt(Sigma),
             lwd = 2, col = cols[i])
}
legend("topleft", col = cols, lty = 1:3, lwd = 2.6, legend = paste0("Sigma = ", Sigmas))