This problem turns out to be well-known in the frequentist literature. In particular, if we use an impropr prior $\Lambda_0=b_0=0$, the posterior scale hyperparameter for the distribution on $\tau$ is

$$\begin{align} b_n&=\frac{1}{2}\left(y^Ty - \mu_n^T\Lambda_n\mu_n\right)\\ &=\frac 12\left(y^Ty -\mu_n^TX^Ty-y^TX\mu_n+ \mu_n^T\Lambda_n\mu_n\right)\\ &=\frac 12\left(y-X\mu_n\right)^T\left(y-X\mu_n\right), \end{align}$$

where we have used the fact that $\Lambda_n=X^TX$ and $\mu_n=\left(X^TX\right)^{-1}X^Ty$. Thus, $b_n$ is $n/2$ times the sample variance of the residuals $y-X\mu_n$. But we already used the data to estimate the regression coefficients and the sample variance of the residuals is a biased estimator of the population variance. In particular, we have $\nu=n-p$ degrees of freedom and an unbiased estimate of the population variance is $$ \frac{1}{n-p}\left(y-X\mu_n\right)^T\left(y-X\mu_n\right). $$

Whenever $n\gg p$ just using the sample variance is fine because $\frac{n}{n-p}\approx 1$. However, as soon as $p$ becomes comparable with $n$ the population variance is underestimated by the sample variance. The inference fails.

I like conjugate priors. But in this case they caused me some trouble. Moral of the story: don't just pick the functional form of your priors because they are conjugate.

As an aside: using a variational mean-field approximation such that the posterior factorises with respect to the regression coefficients and the noise precision works much better than the closed form solution provided by the conjugate priors (better being defined as the posterior being consistent with the true value that was used to generate the data).

This problem turns out to be well-known in the frequentist literature. In particular, if we use an impropr prior $\Lambda_0=b_0=0$, the posterior scale hyperparameter for the distribution on $\tau$ is

$$\begin{align} b_n&=\frac{1}{2}\left(y^Ty - \mu_n^T\Lambda_n\mu_n\right)\\ &=\frac 12\left(y^Ty -\mu_n^TX^Ty-y^TX\mu_n+ \mu_n^T\Lambda_n\mu_n\right)\\ &=\frac 12\left(y-X\mu_n\right)^T\left(y-X\mu_n\right), \end{align}$$

where we have used the fact that $\Lambda_n=X^TX$ and $\mu_n=\left(X^TX\right)^{-1}X^Ty$. Thus, $b_n$ is $n/2$ times the sample variance of the residuals $y-X\mu_n$. But we already used the data to estimate the regression coefficients and the sample variance of the residuals is a biased estimator of the population variance. In particular, we have $\nu=n-p$ degrees of freedom and an unbiased estimate of the population variance is $$ \frac{1}{n-p}\left(y-X\mu_n\right)^T\left(y-X\mu_n\right). $$

Whenever $n\gg p$ just using the sample variance is fine because $\frac{n}{n-p}\approx 1$. However, as soon as $p$ becomes comparable with $n$ the population variance is underestimated by the sample variance. The inference fails. The problem is discussed in the context of maximum likelihood inference on page 388 in "Data Analysis Using Regression and Multilevel/Hierarchical Models".

I like conjugate priors. But in this case they caused me some trouble. Moral of the story: don't just pick the functional form of your priors because they are conjugate.

As an aside: using a variational mean-field approximation such that the posterior factorises with respect to the regression coefficients and the noise precision works much better than the closed form solution provided by the conjugate priors (better being defined as the posterior being consistent with the true value that was used to generate the data).