Normal-Inverse-Gamma Regression With Unknown Variance

Why Model the Variance as Unknown?

In Gaussian regression, treating the error variance as known can make calculations convenient, but in most applications that variance must be learned from the data as well. A conjugate prior for the coefficients and variance preserves tractable posterior distributions: the updated distribution has the same Normal-Inverse-Gamma form as the prior. Its marginal distributions give a Student $t$ distribution for each coefficient and a Student $t$ predictive distribution for a future outcome.

The calculations below are conditional on a specified regression design matrix. They do not resolve whether the chosen regressors or the Gaussian error model are appropriate; those are substantive modeling assumptions to assess before interpreting the posterior.

Regression Model and Prior

Let $y$ be an $n$-vector of outcomes and $X$ an $n\times k$ design matrix, including an intercept if the model requires one. Write $\beta$ for the $k$-vector of regression coefficients and $\sigma^2$ for the error variance. Assume that, conditional on $X,\beta,\sigma^2$, the outcomes have independent normal errors with common variance:

y\mid \beta,\sigma^2,X\sim N(X\beta,\sigma^2 I_n).

Assume $X$ has full column rank. Use the following proper prior: conditional on $\sigma^2$, the coefficient vector is normal, and the variance has an inverse-gamma distribution. The parameters $b_0$, $V_0$, $a_0$, and $d_0$ are fixed, with $V_0$ positive definite and $a_0,d_0>0$.

Definition: Under the Normal-Inverse-Gamma prior, $$ \beta\mid \sigma^2\sim N(b_0,\sigma^2 V_0),\qquad \sigma^2\sim \operatorname{IG}(a_0,d_0). $$ Here $\operatorname{IG}(a,d)$ has density proportional to $(\sigma^2)^{-a-1}\exp(-d/\sigma^2)$ for $\sigma^2>0$; its normalized density is $d^a\{\Gamma(a)\}^{-1}(\sigma^2)^{-a-1}\exp(-d/\sigma^2)$. This parameterization uses $d$ as a scale in the reciprocal-variance exponent.

The prior couples coefficient uncertainty to the error variance: the conditional covariance of $\beta$ is $\sigma^2V_0$. Consequently, learning about $\sigma^2$ affects the marginal uncertainty about $\beta$. This is the key difference from a prior that specifies coefficient uncertainty independently of the variance.

Deriving the Joint Posterior

Define the updated quantities

V_n=(V_0^{-1}+X'X)^{-1},\qquad b_n=V_n(V_0^{-1}b_0+X'y), $$ $$ a_n=a_0+\frac n2,\qquad d_n=d_0+\frac12\left(y'y+b_0'V_0^{-1}b_0-b_n'V_n^{-1}b_n\right).

The posterior is again Normal-Inverse-Gamma, with these updated parameters. The central algebraic step is completing the square in $\beta$. The likelihood times the conditional prior contains the quadratic form

(y-X\beta)'(y-X\beta)+(\beta-b_0)'V_0^{-1}(\beta-b_0).

Expanding in $\beta$, the coefficient of its quadratic term is $V_n^{-1}=X'X+V_0^{-1}$, and its linear term is $-2\beta'(X'y+V_0^{-1}b_0)$. Since $V_n^{-1}b_n=X'y+V_0^{-1}b_0$, completing the square gives

(y-X\beta)'(y-X\beta)+(\beta-b_0)'V_0^{-1}(\beta-b_0) =(\beta-b_n)'V_n^{-1}(\beta-b_n) +y'y+b_0'V_0^{-1}b_0-b_n'V_n^{-1}b_n.

The likelihood contributes a factor $(\sigma^2)^{-n/2}$, and the prior density for $\sigma^2$ contributes $(\sigma^2)^{-a_0-1}$. Their combined power is therefore $(\sigma^2)^{-a_n-1}$; the remaining exponential term has inverse-gamma scale $d_n$. Thus the joint posterior is

Theorem: Under the stated likelihood and proper prior, $$ \beta\mid \sigma^2,y,X\sim N(b_n,\sigma^2V_n),\qquad \sigma^2\mid y,X\sim\operatorname{IG}(a_n,d_n). $$ Equivalently, the joint posterior density is proportional to $$ (\sigma^2)^{-a_n-1-k/2} \exp\left\{-\frac{d_n+\frac12(\beta-b_n)'V_n^{-1}(\beta-b_n)}{\sigma^2}\right\}. $$

The second distribution in the theorem is the marginal posterior for $\sigma^2$, after integrating out $\beta$. Its shape parameter remains $a_n=a_0+n/2$: integrating the conditional normal over its $k$ coefficients cancels the corresponding $k/2$ power in the joint density.

Worked Example: Updating the Joint Posterior

Hypothetical data: suppose three outcomes are $1,2,3$, and the model has only an intercept. Thus $X=(1,1,1)'$, $y=(1,2,3)'$. Set $b_0=0$, $V_0=1$, and $\sigma^2\sim\operatorname{IG}(2,1)$. The errors are assumed independent, normal, and of common variance, and the specified prior is proper.

Here $X'X=3$, $X'y=6$, and $y'y=14$. The updated coefficient variance matrix and mean are

$V_n=(1+3)^{-1}=0.25$, and $b_n=0.25(6)=1.5$.

The variance parameters are $a_n=2+3/2=3.5$ and $d_n=1+\frac12(14-1.5^2/0.25)=1+\frac12(14-9)=3.5$. Thus $\beta\mid\sigma^2,y\sim N(1.5,0.25\sigma^2)$, while $\sigma^2\mid y\sim\operatorname{IG}(3.5,3.5)$. These posterior parameters summarize the joint update; the marginal coefficient distribution follows by integrating over the unknown variance.

Marginal Coefficient Distribution and Credible Intervals

Integrating $\sigma^2$ out of the conditional normal posterior produces a multivariate Student $t$ distribution. With the usual scale-matrix convention for the $t$ distribution, its degrees of freedom are $2a_n$, its location is $b_n$, and its scale matrix is $(d_n/a_n)V_n$:

\beta\mid y,X\sim t_{2a_n}\left(b_n,\frac{d_n}{a_n}V_n\right).

For coefficient $j$, let $e_j$ be the vector with a 1 in position $j$ and zeroes elsewhere. Its marginal posterior is a univariate $t$ distribution with location $b_{n,j}$ and scale $s_j=\sqrt{(d_n/a_n)(V_n)_{jj}}$. In particular, a central $100(1-\alpha)\%$ marginal credible interval is

b_{n,j}\ \pm\ t_{2a_n,\,1-\alpha/2}\sqrt{\frac{d_n}{a_n}(V_n)_{jj}},

where $t_{\nu,q}$ denotes the $q$-quantile of a Student $t$ distribution with $\nu$ degrees of freedom. The interval describes posterior uncertainty about that coefficient, conditional on the model, prior, data, and design. It is marginal for one coefficient; it is not automatically a joint credible region for all coefficients.

Worked Example: A Marginal Credible Interval for the Intercept

Continue with the hypothetical intercept-only regression above. The marginal posterior has $2a_n=7$ degrees of freedom, location $1.5$, and scale $\sqrt{(3.5/3.5)(0.25)}=0.5$. For a 95% interval, the $t_7$ 0.975 quantile is approximately $2.3646$.

The interval is $1.5\pm2.3646(0.5)=1.5\pm1.1823$, which gives approximately $(0.3177,2.6823)$. The coefficient is an intercept, measured in the same units as the outcome. Conditional on this model and prior, the posterior probability that the intercept lies between about 0.318 and 2.682 outcome units is 0.95. This is a Bayesian probability statement about the parameter under the posterior, not a claim that the fixed parameter has a 95% frequentist probability of lying in this realized interval.

Posterior Predictive Distribution

A future observation requires accounting for two sources of uncertainty: uncertainty about the regression mean and the new observation’s error. Let $x_*$ be a specified $k$-vector of regressors for a new case, and suppose its outcome follows the same Gaussian model with a new error independent of the observed errors, conditional on the parameters. Then

y_*\mid\beta,\sigma^2,x_*\sim N(x_*'\beta,\sigma^2).

Conditional on $\sigma^2$ and the observed data, integrating over $\beta\mid\sigma^2,y,X$ gives a normal distribution with mean $x_*'b_n$ and variance $\sigma^2(1+x_*'V_nx_*)$. The 1 is the variance of the new error; the other term reflects uncertainty about the regression mean. Integrating over the inverse-gamma posterior for $\sigma^2$ gives

Result: The posterior predictive distribution for a single future outcome is $$ y_*\mid y,X,x_*\sim t_{2a_n}\left(x_*'b_n,\frac{d_n}{a_n}(1+x_*'V_nx_*)\right), $$ where the second argument is the squared scale of the univariate Student $t$ distribution.

For a vector of future outcomes at a specified design matrix $X_*$, with independent new errors conditional on the parameters, the corresponding multivariate $t$ scale matrix is $(d_n/a_n)(I+X_*V_nX_*')$. Predicting a future observation is not the same as estimating its conditional mean: the predictive distribution includes the new error and is therefore wider.

Worked Example: Predicting a Future Outcome

For the same hypothetical intercept-only regression, the new design value is $x_*=1$. The predictive location is $x_*'b_n=1.5$, and its squared scale is \[ \frac{d_n}{a_n}(1+x_*'V_nx_*)=1(1+0.25)=1.25. \] The predictive scale is therefore $\sqrt{1.25}\approx1.1180$, with 7 degrees of freedom. This calculation assumes the future case follows the same intercept-only model and has a new, independent error with the same variance.

A central 95% predictive interval is $1.5\pm2.3646(1.1180)$, approximately $1.5\pm2.6437$, or $(-1.144,4.144)$ outcome units. Under the model and prior, this interval contains a future outcome with posterior predictive probability 0.95. It is wider than the coefficient interval because it includes both uncertainty in the intercept and the variation of a new outcome around that intercept.

Common Pitfall and Practical Checks

A common mistake is to use only the uncertainty in the fitted conditional mean when constructing a predictive interval. In the formula for a future outcome, that would omit the 1 from $1+x_*'V_nx_*$. The resulting interval would describe uncertainty about the mean $x_*'\beta$, not the outcome $y_*$, and would be too narrow for prediction.

Before using these results, verify that the assumptions match the intended analysis: the regression is linear in the coefficients, the design matrix has full column rank, errors are conditionally independent normal with a common variance, and the stated prior is proper with positive-definite $V_0$ and positive $a_0,d_0$. The posterior calculations do not establish that the model is substantively adequate. Also report the outcome units for coefficient and prediction intervals, and remember that the coefficient interval above is marginal rather than simultaneous.

Key takeaway: With a Normal-Inverse-Gamma prior, the posterior updates through $V_n,b_n,a_n,d_n$. Integrating out the unknown variance gives Student $t$ coefficient marginals, while prediction additionally incorporates the new observation’s error variance.

Check Your Understanding

Try these questions using the stated parameterization.

In the inverse-gamma prior $\operatorname{IG}(a_0,d_0)$, which parameter appears in the exponent as $d_0/\sigma^2$?
Why does the posterior predictive variance include the term 1 in addition to $x_*'V_nx_*$?
For a coefficient with posterior location $b_{n,j}$, scale squared $(d_n/a_n)(V_n)_{jj}$, and $2a_n$ degrees of freedom, write its central 95% marginal credible interval.
Does a marginal 95% credible interval for each coefficient necessarily give a 95% joint credible region for the entire coefficient vector? Explain.
What additional assumption about a future observation is needed to use the stated posterior predictive distribution?