Sparse Priors

1 Normal prior

In what follows, \(\theta\) is a coefficient in the real line and the priors we will present will increasingly induce shrinkage and sparsity. We start with the normal prior: \[ \theta \sim N(0,1), \] with p.d.f. (probability density function) \[ p(\theta) = \frac{1}{\sqrt{2\pi}}\exp\{-0.5\theta^2\} \]

M      = 1000000
thetas = seq(-5,5,length=M)
h      = thetas[2]-thetas[1]
theta0 = 10


par(mfrow=c(1,1))
theta.n = rnorm(M)
den.n   = dnorm(thetas)

hist(theta.n,breaks=seq(min(theta.n),max(theta.n),length=100),xlim=c(-5,5),
     prob=TRUE,xlab=expression(theta),main="Normal")
lines(thetas,den.n,col=2,lwd=2)

2 Student’s \(t\) prior

When \[ \theta \sim t_\nu(0,1) \] its p.d.f is given by \[ p(\theta) = \frac{\Gamma\left(\frac{\nu+1}{2}\right)}{\Gamma\left(\frac{\nu}{2}\right)\sqrt{\pi \nu}}\left(1+\theta^2\right)^{-\frac{\nu+1}{2}}. \] The Student’s \(t\) distribution is one of the distribution derived form a scale-mixture of normals, i.e. \[ \mbox{If} \ \ (\theta|\tau^2) \sim N(0,\tau^2) \ \ \mbox{and} \ \ \tau^2 \sim IG(\nu/2,\nu/2), \ \ \mbox{then} \ \ \theta \sim t_\nu(0,1). \]

nu      = 4
theta.t = rt(M,df=nu)
theta.t = theta.t[abs(theta.t)<theta0]
den.t   = dt(thetas,df=nu)

hist(theta.t,breaks=c(min(theta.t),seq(-5,5,length=100),max(theta.t)),xlim=c(-5,5),
     prob=TRUE,xlab=expression(theta),main="Student's t")
lines(thetas,den.t,col=2,lwd=2)

3 Laplace prior

When \[ \theta \sim \mbox{Laplace}(0,1) \] its p.d.f is given by \[ p(\theta) = \frac12\exp\{-|\theta|\}. \] The Laplace distribution, also known as the double-exponential distribution, is another distribution derived form a scale-mixture of normals, i.e. \[ \mbox{If} \ \ (\theta|\tau^2) \sim N(0,\tau^2) \ \ \mbox{and} \ \ \tau^2 \sim G(1,1/2), \ \ \mbox{then} \ \ \theta \sim \mbox{Laplace}(0,1). \]

tau2    = rgamma(M,1,1/2)
theta.l = rnorm(M,0,sqrt(tau2))
den.l   = 0.5*exp(-abs(thetas))

hist(theta.l,breaks=c(min(theta.l),seq(-5,5,length=100),max(theta.l)),xlim=c(-5,5),
     prob=TRUE,xlab=expression(theta),main="Laplace")
lines(thetas,den.l,col=2,lwd=2)

4 Horseshoe prior

The parameter \(\theta\) follows a horseshoe prior when it is derived according to the following scheme: \[ \theta | \tau^2 \sim N(0,\tau^2), \ \ \tau|\lambda \sim C^{+}(0,\lambda) \ \ \mbox{and} \ \ \lambda \sim C^{+}(0,1), \] where \(C^{+}(0,\lambda)\) is the half-Cauchy distribution with p.d.f \[ p(\tau|\lambda) = \frac{2}{\pi}\left(1+\frac{\tau^2}{\lambda^2}\right)^{-1}. \] The horseshoe density can be approximated by \[ p(\theta) \approx \log\left(1+\frac{2}{\theta^2}\right). \]

lambda  = abs(rt(M,1))
tau     = abs(sqrt(lambda)*rt(M,1))
theta.h = rnorm(M,0,tau)
theta.h = theta.h[abs(theta.h)<theta0]

hist(theta.h,breaks=c(min(theta.h),seq(-5,5,length=100),max(theta.h)),xlim=c(-5,5),
     prob=TRUE,xlab=expression(theta),main="Horseshoe")
den.h = log(1+2/thetas^2)
k = sum(den.h*h)
den.h = den.h/k
lines(thetas,den.h,col=2,lwd=2)

5 Normal-Gamma prior

When \[ \theta \sim \mbox{Normal-Gamma}(\lambda,\gamma) \] its p.d.f is given by \[ p(\theta) = \frac{1}{\sqrt{\pi}2^{\lambda-1/2}\gamma^{\lambda+1/2}\Gamma(\lambda)}|\theta|^{\lambda-1/2}K_{\lambda-1/2}\left(\frac{|\theta|}{\gamma}\right), \] onde \(K\) is the modified Bessel function of the third kind.

The Normal-Gamma distribution is another distribution derived form a scale-mixture of normals, i.e. \[ \mbox{If} \ \ (\theta|\tau^2) \sim N(0,\tau^2) \ \ \mbox{and} \ \ \tau^2 \sim G(\lambda,1/(2\gamma^2)), \ \ \mbox{then} \ \ \theta \sim \mbox{Normal-Gamma}(\lambda,\gamma). \] It can be shown that \(V(\theta) = 2\lambda \gamma^2\) and the excess kurtosis is equal to \(3/\lambda\), so \(\lambda\) controls the heaviness of the tails. The Laplace\((0,1)\) is the special case when \(\lambda=\gamma=1\). Below we set \(\gamma=1\) and \(\lambda=0.5\), such that \(V(\theta)=1\) and the excess kurtosis is equal to 6.

lambda   = 0.5
gamma    = 1
tau2     = rgamma(M,lambda,1/(2*gamma^2))
theta.ng = rnorm(M,0,sqrt(tau2))

hist(theta.ng,breaks=c(min(theta.ng),seq(-5,5,length=100),max(theta.ng)),xlim=c(-5,5),
     prob=TRUE,xlab=expression(theta),main="Normal-Gamma")
const = sqrt(pi)*2^(lambda-0.5)*gamma^(lambda+0.5)*gamma(lambda)
den.ng = (1/const)*abs(thetas)^(lambda-0.5)*besselK(abs(thetas)/gamma,lambda-0.5)
lines(thetas,den.ng,col=2,lwd=2)

6 All together

par(mfrow=c(1,2))

plot(thetas,den.t,xlim=c(-5,5),ylim=c(0,2),lwd=1.5,xlab=expression(theta),main="",type="l",ylab="Density")
lines(thetas,den.l,col=2,lwd=1.5)
lines(thetas,den.h,col=3,lwd=1.5)
lines(thetas,den.ng,col=4,lwd=1.5)
lines(thetas,dnorm(thetas,0,1),col=6,lwd=1.5)
legend("topright",legend=c("Normal","Student's t","Laplace","Horseshoe","Normal-Gamma"),col=c(6,1:4),
lty=1,lwd=1.5,bty="n")

plot(thetas,log(den.t),xlim=c(-5,5),ylim=c(-7,2),lwd=1.5,xlab=expression(theta),main="",type="l",ylab="Log density")
lines(thetas,log(den.l),col=2,lwd=1.5)
lines(thetas,log(den.h),col=3,lwd=1.5)
lines(thetas,log(den.ng),col=4,lwd=1.5)
lines(thetas,dnorm(thetas,0,1,log=TRUE),col=6,lwd=1.5)
legend("topright",legend=c("Normal","Student's t","Laplace","Horseshoe","Normal-Gamma"),col=c(6,1:4),
lty=1,lwd=1.5,bty="n")