# 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. – jbowman Dec 17 '12 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. – far away Dec 17 '12 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.) – jbowman Dec 17 '12 at 19:55
• Thank you for your comments. I edited the question according to them. I hope that it makes a better sense now. – far away Dec 17 '12 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. – jbowman Dec 17 '12 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.