Bayesian linear regression

Consider standard linear regression problem, where we specify the conditional density of y given x predictor variables:

$y_{i}=\beta x_{i}+\epsilon _{i},\,$

where the noise $\epsilon$ is i.i.d. and normally distributed

$\epsilon _{i}\sim N(0,\sigma ^{2}).\,$

A common, least linear squares solution, is to estimate the slope ${\hat {\beta }}$ using the Moore-Penrose pseudoinverse:

${\hat {\beta }}=(X^{T}X)^{-1}X^{T}y$.

where $X$ is the vector of $x_{i}$ (of length $n$).

This is a frequentists view, and assumes we have enough measurements of $x_{i}$ to say something meaningful about y. In the Empirical Bayes approach, we will assume we have only a small sample of $x_{i}$ for our individual measurement, and we seek to correct our estimate by "borrowing" information from a larger set of similar observations.

Let us write our conditional likelihood as

$\rho (y|X,\beta ,\sigma ^{2})\propto (\sigma ^{2})^{-n/2}exp(-{\frac {1}{2{\sigma }^{2}}}(y-\beta X)^{T}(y-\beta X)),\,$

We seek a natural conjugate prior--a joint density $\rho (\beta ,\sigma ^{2})$ which is of the same functional form as the likelihood. Since the likelihood is quadratic in $\beta$, we re-write the likelihood so it is normal in $(\beta -{\hat {\beta }})$. Write

$(y-\beta X)^{T}(y-\beta X)=(y-{\hat {\beta }}X)^{T}(y-{\hat {\beta }}X)+(\beta -{\hat {\beta }})^{T}(X^{T}X)(\beta -{\hat {\beta }})$

Now re-write the likelihood as

$\rho (y|X,\beta ,\sigma ^{2})\propto (\sigma ^{2})^{-v/2}exp(-{\frac {vs^{2}}{2{\sigma }^{2}}})(\sigma ^{2})^{-(n-v)/2}exp(-{\frac {1}{2{\sigma }^{2}}}(\beta -{\hat {\beta }})^{T}(X^{T}X)(\beta -{\hat {\beta }})),\,$

where

$vs^{2}=(y-{\hat {\beta }}X)^{T}(y-{\hat {\beta }}X)\ ,v=n-k$

This suggests a form for the priors:

$\rho (\beta ,\sigma ^{2})=\rho (\sigma ^{2})\rho (\beta |\sigma ^{2}),\,$

where $\rho (\sigma ^{2})$ is an inverse-gamma distribution

$\rho (\sigma ^{2})\propto (\sigma ^{2})^{-v_{0}/2+1}exp(-{\frac {v_{0}s_{0}^{2}}{2{\sigma }^{2}}}),\,$

and $\rho (\beta |\sigma ^{2})$ is an Normal distribution

$\rho (\beta |\sigma ^{2})\propto (\sigma ^{2})^{-k}exp(-{\frac {1}{2{\sigma }^{2}}}(\beta -{\bar {\beta }})^{T}(A)(\beta -{\bar {\beta }})),\,$

With the prior now specified, we can express the posterior distribution as

$\rho (\beta ,\sigma ^{2}|y,X)\propto \rho (y|X,\beta ,\sigma ^{2})\rho (\beta |\sigma ^{2})\rho (\sigma ^{2})$

$\propto (\sigma ^{2})^{-n/2}exp(-{\frac {1}{2{\sigma }^{2}}}(y-\beta X)^{T}(y-\beta X))$

$\times (\sigma ^{2})^{-k}exp(-{\frac {1}{2{\sigma }^{2}}}(\beta -{\bar {\beta }})^{T}(A)(\beta -{\bar {\beta }})).$

$\times (\sigma ^{2})^{-v_{0}/2+1}exp(-{\frac {v_{0}s_{0}^{2}}{2{\sigma }^{2}}})$

With some re-arrangement, we can re-write the posterior so that the posterior mean ${\tilde {\beta }}$ is weighted average of the least squares estimator and the prior mean:

${\tilde {\beta }}=(X^{T}X+A)^{-1}(X^{T}X{\hat {\beta }}+A{\bar {\beta }})$

where $U$ comes from the LU decomposition of $A$ (which is a positive-definite matrix by design)

$A=U^{T}U.$

This is the key result of the Empirical Bayes approach; it allows us to estimate the slope $\beta$ for our original linear regression problem by combining estimates using the least squares estimate ${\hat {\beta }}$ for a single set of measurements with the empirical prior estimate ${\bar {\beta }}$ from a large collection of similar measurements. (Notice that the weighted average also depends on the empirical estimate of the prior covariance matrix $A$.)

To justify this, collect the quadratic terms in the exponential and and now express this as a quadratic form in $\beta -{\tilde {\beta }}$:

$(y-\beta X)^{T}(y-\beta X))+(\beta -{\bar {\beta }})^{T}(A)(\beta -{\bar {\beta }})=(v-W\beta )^{T}(v-W\beta )$

$=ns^{2}+(\beta -{\bar {\beta }})^{T}W^{T}W(\beta -{\bar {\beta }})$

where

$ns^{2}=(v-W{\tilde {\beta }})^{T}(v-W{\tilde {\beta }}),v=[y,U{\bar {B}}],W=[X,U]$

The posterior can now be expressed as a Normal distribution $N({\tilde {\beta }},\sigma ^{2}(X^{T}X+A)^{-1}$ times an inverse-gamma distribution:

$\rho (\beta ,\sigma ^{2}|y,X)\propto (\sigma ^{2})^{-k/2}exp(-{\frac {1}{2{\sigma }^{2}}}(\beta -{\bar {\beta }})^{T}(X^{T}X+A)(\beta -{\bar {\beta }}))\times (\sigma ^{2})^{-(n+v_{0})/2+1}exp(-{\frac {(v_{0}s_{0}^{2}+ns^{2})}{2{\sigma }^{2}}})$

A similar analysis can be performed for general case of multi-variate regression for a Bayesian Estimation of covariance matrices.

Example:

Suppose the weights of a large population of 35-year-old men are normally distributed with expected value μ and standard deviation σ. A crude measuring instrument measures a man's weight with a measurement error that is normally distributed with expected value 0 and standard deviation τ. The man's true weight is not observable; his weight measured with error is observed. The conditional probability distribution of a randomly chosen man's true weight, given his weight-measured-with-error, can be found by using Bayes' theorem, and then the conditional expected value can be used as an estimate of his true weight, provided that the values of μ, σ, and τ are known. But they are not. One may use the data to estimate the standard deviation of the measurement errors by measuring each man multiple times. One may similarly estimate the population average weight and the population standard deviation of weights by weighing multiple men. These estimates of parameters based on the data are the occasion for the use of the word empirical. Finally, one may then estimate the aforementioned conditional expected true weight by using Bayes' theorem.

References

• Bradley P. Carlin and Thomas A. Louis, Bayes and Empirical Bayes Methods for Data Analysis, Chapman & Hall/CRC, Second edition 2000,

• Peter E. Rossi, Greg M. Allenby, and Robert McCulloch, Bayesian Statistics and Marketing, John Wiley & Sons, Ltd, 2006