# Bayes regression: how is it done in comparison to standard regression?

I got some questions about the Bayesian regression:

1. Given a standard regression as $y = \beta_0 + \beta_1 x + \varepsilon$. If I want to change this into a Bayesian regression, do I need prior distributions both for $\beta_0$ and $\beta_1$ (or doesn't it work this way)?

2. In standard regression one would try to minimize the residuals to get single values for $\beta_0$ and $\beta_1$. How is this done in Bayes regression?

I really struggle a lot here:

$$\text{posterior} = \text{prior} \times \text{likelihood}$$

Likelihood comes from the current dataset (so it's my regression parameter but not as a single value but as a likelihood distribution, right?). Prior comes from a previous research (let's say). So I got this equation:

$$y = \beta_1 x + \varepsilon$$

with $\beta_1$ being my likelihood or posterior (or is this just totally wrong)?

I simply can't understand how the standard regression transforms into a Bayes one.

## 2 Answers

Simple linear regression model

$$y_i = \alpha + \beta x_i + \varepsilon$$

can be written in terms of probabilistic model behind it

$$\mu_i = \alpha + \beta x_i \\ y_i \sim \mathcal{N}(\mu_i, \sigma)$$

i.e. dependent variable $Y$ follows normal distribution parametrized by mean $\mu_i$, that is a linear function of $X$ parametrized by $\alpha,\beta$, and by standard deviation $\sigma$. If you estimate such model using ordinary least squares, you do not have to bother about the probabilistic formulation, because you are searching for optimal values of $\alpha,\beta$ parameters by minimizing the squared errors of fitted values to predicted values. On another hand, you could estimate such model using maximum likelihood estimation, where you would be looking for optimal values of parameters by maximizing the likelihood function

$$\DeclareMathOperator*{\argmax}{arg\,max} \argmax_{\alpha,\,\beta,\,\sigma} \prod_{i=1}^n \mathcal{N}(y_i; \alpha + \beta x_i, \sigma)$$

where $\mathcal{N}$ is a density function of normal distribution evaluated at $y_i$ points, parametrized by means $\alpha + \beta x_i$ and standard deviation $\sigma$.

In Bayesian approach instead of maximizing the likelihood function alone, we would assume prior distributions for the parameters and use Bayes theorem

$$\text{posterior} \propto \text{likelihood} \times \text{prior}$$

The likelihood function is the same as above, but what changes is that you assume some prior distributions for the estimated parameters $\alpha,\beta,\sigma$ and include them into the equation

$$\underbrace{f(\alpha,\beta,\sigma\mid Y,X)}_{\text{posterior}} \propto \underbrace{\prod_{i=1}^n \mathcal{N}(y_i\mid \alpha + \beta x_i, \sigma)}_{\text{likelihood}} \; \underbrace{f_{\alpha}(\alpha) \, f_{\beta}(\beta) \, f_{\sigma}(\sigma)}_{\text{priors}}$$

"What distributions?" is a different question, since there is unlimited number of choices. For $\alpha,\beta$ parameters you could, for example assume normal distributions parametrized by some hyperparameters, or $t$-distribution if you want to assume heavier tails, or uniform distribution if you do not want to make much assumptions, but you want to assume that the parameters can be a priori "anything in the given range", etc. For $\sigma$ you need to assume some prior distribution that is bounded to be greater then zero, since standard deviation needs to be positive. This may lead to the model formulation as illustrated below by John K. Kruschke. While in maximum likelihood you were looking for a single optimal value for each of the parameters, in Bayesian approach by applying Bayes theorem you obtain the posterior distribution of the parameters. The final estimate will depend on the information that comes from your data and from your priors, but the more information is contained in your data, the less influential are priors.

Notice that when using uniform priors, they take form $f(\theta) \propto 1$ after dropping the normalizing constants. This makes Bayes theorem proportional to likelihood function alone, so the posterior distribution will reach it's maximum at exactly the same point as maximum likelihood estimate. What follows, the estimate under uniform priors will be the same as by using ordinary least squares since minimizing the squared errors corresponds to maximizing the normal likelihood.

To estimate a model in Bayesian approach in some cases you can use conjugate priors, so the posterior distribution is directly available (see example here). However in vast majority of cases posterior distribution will not be directly available and you will have to use Markov Chain Monte Carlo methods for estimating the model (check this example of using Metropolis-Hastings algorithm to estimate parameters of linear regression). Finally, if you are only interested in point estimates of parameters, you could use maximum a posteriori estimation, i.e.

$$\argmax_{\alpha,\,\beta,\,\sigma} f(\alpha,\beta,\sigma\mid Y,X)$$

For more detailed description of logistic regression you can check the Bayesian logit model - intuitive explanation? thread.

For learning more you could check the following books:

Kruschke, J. (2014). Doing Bayesian Data Analysis: A Tutorial with R, JAGS, and Stan. Academic Press.

Gelman, A., Carlin, J. B., Stern, H. S., and Rubin, D. B. (2004). Bayesian data analysis. Chapman & Hall/CRC.

• +1 Given the way the question is stated, I'd maybe emphasize a bit more this philosophical difference: In ordinary least squares and maximum likelihood estimation, we are starting with the question "What are the best values for $\beta_i$ (perhaps for later use)?", whereas in the full Bayesian approach, we start with the question "What can we say about the unknown values $\beta_i$?" and then maybe proceed to using the maximum a posteriori or posterior mean if a point estimate is needed. – JiK Dec 22 '16 at 11:32
• +1. One more thing that might be useful to point out to clarify the relationship between Bayesian and OLS approaches is that OLS can be understood as posterior mean under a flat prior (at least as far as I understand). Would be great if you could elaborate on that a bit in your answer. – amoeba Dec 22 '16 at 15:07
• @amoeba it's a good point, I'll think about it. But on another hand, I don't want to make the answer overtly long, so there is a point in going to details. – Tim Dec 22 '16 at 16:08
• @amoeba FYI, I added a brief comment on that. – Tim Dec 28 '16 at 11:52

Given a data set $D = (x_1,y_1), \ldots, (x_N,y_N)$ where $x \in \mathbb{R}^d, y \in \mathbb{R}$, a Bayesian Linear Regression models the problem in the following way:

Prior: $$w \sim \mathcal{N}(0, \sigma_w^2 I_d)$$

$w$ is vector $(w_1, \ldots, w_d)^T$, so the previous distribution is a multivariate Gaussian; and $I_d$ is the $d\times d$ identity matrix.

Likelihood: $$Y_i \sim \mathcal{N}(w^T x_i, \sigma^2)$$

We assume that $Y_i \perp Y_j | w, i \neq j$

For now we'll use the precision instead of the variance, $a = 1/\sigma^2$, and $b = 1/\sigma_w^2$. We'll also assume that $a,b$ are known.

The prior can be stated as $$p(w) \propto \exp \Big\{ -\frac{b}{2} w^t w \Big\}$$

And the likelihood $$p(D|w) \propto \exp \Big\{ -\frac{a}{2} (y-Aw)^T (y-Aw) \Big\}$$

where $y = (y_1,\ldots,y_N)^T$ and $A$ is a $n\times d$ matrix where the i-th row is $x_i^T$.

Then the posterior is $$p(w|D) \propto p(D|w) p(w)$$

After many calculations we discover that

$$p(w|D) \sim \mathcal{N}(w | \mu, \Lambda^{-1})$$

where ($\Lambda$ is the precision matrix)

$$\Lambda = a A^T A + b I_d$$ $$\mu = a \Lambda^{-1} A^T y$$

Notice that $\mu$ is equal to the $w_{MAP}$ of the regular linear regression, this is because for the Gaussian, the mean is equal to the mode.

Also, we can make some algebra over $\mu$ and get the following equality ($\Lambda = aA^TA+bI_d$):

$$\mu = (A^T A + \frac{b}{a} I_d)^{-1} A^T y$$

and compare with $w_{MLE}$:

$$w_{MLE} = (A^T A)^{-1} A^T y$$

The extra expression in $\mu$ corresponds to the prior. This is similar to the expression for the Ridge regression, for the special case when $\lambda = \frac{b}{a}$. Ridge regression is more general because the technique can choose improper priors (in the Bayesian perspective).

For the predictive posterior distribution:

$$p(y|x,D) = \int p(y|x,D,w) p(w|x,D) dw = \int p(y|x,w) p(w|D) dw$$

it is possible to calculate that

$$y|x,D \sim \mathcal{N}(\mu^Tx, \frac{1}{a} + x^T \Lambda^{-1}x)$$

Reference: Lunn et al. The BUGS Book

For using a MCMC tool like JAGS/Stan check Kruschke's Doing Bayesian Data Analysis

• Thank you jpneto. i feel that this is a great answer, but i dont understand it yet because of a lack in math-knowledge. But i will definitely read it again after gaining some math-skills – TinglTanglBob Dec 21 '16 at 18:53
• This is very nice, but the assumption that the precision is known is a bit uncommon. Isn't it much more common to assume an inverse Gamma distribution for the variance, i.e., a Gamma distribution for the precision? – DeltaIV Dec 22 '16 at 10:57
• +1. Can you comment a bit more on "Ridge regression is more general because the technique can choose improper priors"? I don't get it. I thought RR = Gaussian (proper) prior on $w$. – amoeba Dec 23 '16 at 8:05
• @amoeba: The Gaussian prior is $w \sim N(0,\lambda^{-1} I_d)$ but $\lambda$ can be zero which results on an improper prior, ie, it results in the MLE. – jpneto Dec 23 '16 at 9:25
• @DeltaIV: sure, when we have uncertainty about a parameter we can model that with a prior. The assumption of known precision is to make it easier to find an analytic solution. Usually, those analytic solutions are not possible and we must use approximations, like MCMC or some variational technique. – jpneto Dec 23 '16 at 9:29