# Laplace distribution
https://statisticaloddsandends.wordpress.com/2018/12/21/laplace-distribution-as-a-mixture-of-normals/

t4 = function(theta){
  sigma = sqrt(1/2)
  dt(theta/sigma,df=4)/sigma
}

laplace = function(theta){
  exp(-abs(theta))/2
}

horseshoe = function(theta,k){
  log(1+2/theta^2)/k
}

normalgamma = function(theta,lambda){
  gamma = sqrt(1/(2*lambda))	
  f = 1/(sqrt(pi)*2^(lambda-1/2)*(gamma)^(lambda+1/2)*gamma(lambda))*(abs(theta))^(lambda-1/2)*besselK(abs(theta)/gamma, nu=lambda-1/2)
}


theta = seq(-15,15,length=1000)

k = sum(horseshoe(theta,1)*(theta[2]-theta[1]))


par(mfrow=c(1,2))
plot(theta,horseshoe(theta,k),type="l",lwd=2,ylim=c(0,1),xlim=c(-5,5),xlab=expression(theta),ylab="Density")
lines(theta,normalgamma(theta,0.1),col=2,lwd=2)
lines(theta,laplace(theta),col=3,lwd=2)
lines(theta,t4(theta),col=4,lwd=2)
lines(theta,dnorm(theta),col=5,lwd=2)
legend("topleft",legend=c("N(0,1)","t(0,1/4,4)","Laplace","Normal-Gamma","Horseshoe"),col=5:1,lwd=2,bty="n",lty=1)

plot(theta,log(horseshoe(theta,k)),type="l",ylim=c(-20,0),lwd=2,xlab=expression(theta),ylab="Log density")
lines(theta,log(normalgamma(theta,0.1)),col=2,lwd=2)
lines(theta,log(laplace(theta)),col=3,lwd=2)
lines(theta,log(t4(theta)),col=4,lwd=2)
lines(theta,dnorm(theta,log=TRUE),col=5,lwd=2)




set.seed(31415)
M = 1000000
par(mfrow=c(1,3))

# Laplace
w = rgamma(M,1,1/2)
theta.l = rnorm(M,0,sqrt(w))

breaks=seq(min(theta.l),max(theta.l),length=100)
hist(theta.l,prob=TRUE,xlab=expression(theta),xlim=c(-10,10),breaks=breaks,main="Laplace draws")
legend("topleft",legend=paste("% (-10,10) : ",round(mean(abs(theta.l)<10),6),sep=""),bty="n")

# Horseshoe
tau    = abs(rt(M,df=1))
lambda = abs(tau*rt(M,df=1))
theta.h = rnorm(M,0,lambda)
breaks=c(min(theta.h),seq(-10,10,length=100),max(theta.h))
hist(theta.h,prob=TRUE,xlab=expression(theta),xlim=c(-10,10),breaks=breaks,main="Horseshoe draws")
legend("topleft",legend=paste("% (-10,10) : ",round(mean(abs(theta.h)<10),6),sep=""),bty="n")



# Normal-Gamma
lambda   = 0.1
gamma    = sqrt(1/(2*lambda))
tau2     = rgamma(M,lambda,1/(2*gamma^2))
theta.ng = rnorm(M,0,sqrt(tau2))

breaks=seq(min(theta.ng),max(theta.ng),length=100)
hist(theta.ng,prob=TRUE,xlab=expression(theta),xlim=c(-10,10),breaks=breaks,main="Normal-Gamma draws")
legend("topleft",legend=paste("% (-10,10) : ",round(mean(abs(theta.ng)<10),6),sep=""),bty="n")





# Normal-Gamma
lambdas = c(0.01,0.05,0.1,0.5,0.75,1)
par(mfrow=c(2,3))
for (lambda in lambdas){
  gamma    = sqrt(1/(2*lambda))
  tau2     = rgamma(M,lambda,1/(2*gamma^2))
  theta.ng = rnorm(M,0,sqrt(tau2))
  breaks=seq(min(theta.ng),max(theta.ng),length=200)
  hist(theta.ng,prob=TRUE,xlab=expression(theta),xlim=c(-10,10),breaks=breaks,main="")
  legend("topleft",legend=paste("% (-10,10) : ",round(mean(abs(theta.ng)<10),6),sep=""),bty="n")
}