# Computing posterior distribution of bayesian lasso

I have a model:

$$\mathbf{y} = \mu\mathbf{1}_n + \mathbf{X}\boldsymbol{\beta} + \boldsymbol{\epsilon}$$

where $\boldsymbol{\epsilon} \sim N(\mathbf{0},\sigma^2\mathbf{I}_n)$.

I have a joint prior: $$\pi(\boldsymbol{\beta}, \sigma^2, \mu) = \pi(\mu) \pi(\sigma^2) \prod\limits_{j=1}^{p}\frac{\lambda}{2\sqrt{\sigma^2}}e^{-\lambda|\beta_j|/\sqrt{\sigma^2}}$$

where $\pi(\mu) \propto 1$.

I want to compute joint posterior $\pi(\boldsymbol{\beta}, \sigma^2, \mu |\mathbf{y})$ and then marginalize out $\mu$.

MY SOLUTION:

According to me, likelihood function is: $$f(\mathbf{y} |\boldsymbol{\beta}, \sigma^2, \mu) = \frac{1}{(2\pi\sigma^2)^{\frac{n}{2}}} \exp\left( -\frac{1}{2\sigma^2}\parallel \mathbf{y}- \mu \mathbf{1}_n -\mathbf{X}\boldsymbol{\beta}\parallel_2^2 \right)$$

Using the Bayes theorem, I got the conclusion that joint log-posterior is proportional to: $$\underbrace{\ln[\pi(\mu)]}_{\to \text{const.}} + \ln[\pi(\sigma^2)]-\frac{p+n}{2}\ln(\sigma^2) - \lambda \parallel \beta \parallel_1 \frac{1}{\sqrt{\sigma^2}} - \frac{1}{2\sigma^2}\parallel \mathbf{y}- \mu\mathbf{1}_n - \mathbf{X}\boldsymbol{\beta} \parallel_2^2$$

HOWEVER, in paper about Bayessian Lasso (2008) written by Park and Casella they got: $$\ln[\pi(\sigma^2)]-\frac{p+n-1}{2}\ln(\sigma^2) - \lambda \parallel \beta \parallel_1 \frac{1}{\sqrt{\sigma^2}} - \frac{1}{2\sigma^2}\parallel \mathbf{\tilde{y}}-\mathbf{X}\boldsymbol{\beta} \parallel_2^2$$

where $\mathbf{\tilde{y}} = \mathbf{y} - \bar{y}\mathbf{1}_n$

Can someone tell me, what is the next step which will lead to a nice marginalization of $\mu$ and getting their result?

• Typically, $X$ will have a constant term in it, and you might not want to shrink the coefficient of that towards 0... think about what that would imply for your prior. Dec 17, 2012 at 18:02
• OK, then lets say that $\mathbf{y}$ is demeaned by overall mean $\bar{y}$ and $\mathbf{X}$ does not include constant term. Dec 17, 2012 at 18:43
• If you do that, then your posterior is conditional upon $\mu = \bar{y}$, not what you want. You want to integrate out $\mu$ instead. Note that in the original paper it starts out with $\mu$ (separated from $X$), but on p. 683, column 2, they refer to "Marginalizing over $\mu$..." and the $n-1$ appears shortly thereafter. For an analogous take on it, in classical statistics, how many degrees of freedom do you have left after you've subtracted off the sample mean? (Think of t tests.) Dec 17, 2012 at 19:55
• Thank you for your comments. I edited the question according to them. I hope that it makes a better sense now. Dec 17, 2012 at 20:59
• Excellent! Now, when you're coming up with the joint log-posterior, it has $\mu$ in it, but Park and Casella's version doesn't, because it's the log of the marginal posterior (with $\mu$ integrated out.) So you'll have to integrate $\mu$ out of your expression for $f(y|\beta, \sigma^2, \mu)$, then take the log to get the "log marginal posterior". You may (or may not) know how to do this; more hints are forthcoming if you don't. Dec 17, 2012 at 21:13

So, we have:

1) Joint prior: $$\pi(\mu, \sigma^2, \boldsymbol{\beta}) = \pi(\sigma^2)\pi(\mu)\prod\limits_{j=1}^{p}\frac{\lambda}{2\sqrt{\sigma^2}}e^{-\lambda|\beta_j|/\sqrt{\sigma^2}},$$

where $\pi(\mu)\propto 1$.

2) Likelihood function: $$f(\mathbf{y}|\mu, \sigma^2, \boldsymbol{\beta}) = \frac{1}{(2\pi\sigma^2)^{\frac{n}{2}}} \exp\left( -\frac{1}{2\sigma^2}\| \mathbf{y}- \mu\mathbf{1}_n - \mathbf{X}\boldsymbol{\beta}\|_2^2 \right)$$

Using Bayes theorem we get that posterior distribution is proportional to: $$\pi(\sigma^2, \boldsymbol{\beta} | \mathbf{\tilde{y}}) \propto \pi(\sigma^2) \exp\left( -\frac{1}{2\sigma^2}\| \mathbf{y}- \mu\mathbf{1}_n - \mathbf{X}\boldsymbol{\beta}\|_2^2 \right) \prod\limits_{j=1}^{p}\frac{\lambda}{2\sqrt{\sigma^2}}e^{-\lambda|\beta_j|/\sqrt{\sigma^2}}$$

Now we want to integrate out $\mu$. For that we need to solve: $$\begin{split} \int_{-\infty}^{\infty} \exp\left( -\frac{1}{2\sigma^2}\| \mathbf{y}- \mu\mathbf{1}_n - \mathbf{X}\boldsymbol{\beta}\|_2^2 \right) d\mu =\\ =\int_{-\infty}^{\infty} \exp\left( -\frac{1}{2\sigma^2} \sum\limits_{i=1}^n\left(\bar{y} -\mu + y_i - \bar{y} - \mathbf{X}_i'\boldsymbol{\beta}\right)^2 \right) d\mu = \\ = \int_{-\infty}^{\infty} \exp\left( -\frac{1}{2\sigma^2} \left[n\left(\bar{y} -\mu \right)^2 + \sum\limits_{i=1}^n\left(y_i - \bar{y} - \mathbf{X}_i'\boldsymbol{\beta}\right)^2 + \right. \right. \\ \left. \left. +2\left(\bar{y} -\mu \right)\sum\limits_{i=1}^n\left(y_i - \bar{y} - \mathbf{X}_i'\boldsymbol{\beta}\right) \right]\right) d\mu \end{split}$$

where $\mathbf{X}_i$ is the $i$-th row of matrix $\mathbf{X}$.

There is a forgotten assumption in the setup that elements of $\mathbf{X}$ are centered. Thanks to this assumption, $\sum\limits_{i=1}^n\left(y_i - \bar{y} - \mathbf{X}_i'\boldsymbol{\beta}\right) = 0$. Because:

\begin{gather*} \sum\limits_{i=1}^n\left(y_i - \bar{y} - \mathbf{X}_i'\boldsymbol{\beta}\right) = \sum\limits_{i=1}^n y_i - n\bar{y} -\sum\limits_{i=1}^n\sum\limits_{j=1}^p (x_{ij}^* - \bar{x}_{\cdot j})\beta_j = \\ = n\bar{y} - n\bar{y} -\sum\limits_{i=1}^n\sum\limits_{j=1}^p (x_{ij}^* - \bar{x}_{\cdot j})\beta_j \end{gather*}

For each $j \in \{1, 2, \dots, p\}$ holds: $$\sum\limits_{i=1}^n (x_{ij}^* - \bar{x}_{\cdot j})\beta_{j} = n\bar{x}_{\cdot j}\beta_{j}-n\bar{x}_{\cdot j}\beta_{j} = 0$$

Now back to our integral. Now we can say that it is equal to: \begin{equation*} \begin{split} \exp\left(-\frac{1}{2\sigma^2}\sum\limits_{i=1}^n\left(y_i - \bar{y} - \mathbf{X}_i'\boldsymbol{\beta}\right)^2\right) \int_{-\infty}^{\infty} \exp\left( -\frac{n}{2\sigma^2} \left(\bar{y} -\mu \right)^2 \right) d\mu \end{split} \end{equation*}

From normal distribution, we know that: $$\int_{-\infty}^{\infty} \frac{\sqrt{n}}{\sqrt{2\pi\sigma^2}} \exp\left( -\frac{n}{2\sigma^2} \left(\bar{y} -\mu \right)^2 \right) d\mu = 1 \ \ \ \Rightarrow$$

$$\int_{-\infty}^{\infty} \exp\left( -\frac{n}{2\sigma^2} \left(\bar{y} -\mu \right)^2 \right) d\mu = \frac{\sqrt{2\pi\sigma^2}}{\sqrt{n}}$$

The final result of posterior distribution after integrating out $\mu$ is: $$\pi(\sigma^2, \boldsymbol{\beta} | \mathbf{\tilde{y}}) \propto \pi(\sigma^2) \exp\left(-\frac{1}{2\sigma^2}\sum\limits_{i=1}^n\left(y_i - \bar{y} - \mathbf{X}_i'\boldsymbol{\beta}\right)^2\right) \frac{\sqrt{2\pi\sigma^2}}{\sqrt{n}} \prod\limits_{j=1}^{p}\frac{\lambda}{2\sqrt{\sigma^2}}e^{-\lambda|\beta_j|/\sqrt{\sigma^2}}$$

By taking log and leaving all terms without $\sigma$ we get the result.