Simulating Horseshoe Distribution

Question

Can someone check whether I correctly simulate the Horseshoe pdf?

lambda = rcauchy(1000, location = 0, scale = 1)

horseshoe = rnorm(1000, 0,lambda)

hist(horseshoe, breaks = 'FD')

Answer 1

The horseshoe distribution is a zero-mean Cauchy-Normal scale mixture of normals (as @jbowman points out). It has one parameter, scale $\tau$; the bigger $\tau$ is, the heavier the tails are (the higher the probability of observing large positive or negative values). 1

Your simulation assumes $\tau = 1$. You also have to fold the Cauchy (by taking the absolute value) as it doesn't make sense to sample from a Normal with negative standard deviation. That's why your code throws warnings and you end up with a sample that's about half NaN's.

The horseshoe distribution has two equivalent specifications: $$ \lambda \sim \text{Half-Cauchy}(0, 1) \\ X \sim \text{Normal}(0, \lambda \cdot \tau) $$ and $$ \lambda \sim \text{Half-Cauchy}(0, \tau) \\ X \sim \text{Normal}(0, \lambda) $$

It's straightforward to use either to generate a random sample. In the plot I also overlay the histogram with the probability density function.

1 C. M. Carvalho, N. G. Polson, and J. G. Scott (2010). The horseshoe estimator for sparse signals. Biometrika, 97(2):465–480

tau <- 2 # scale parameter of horseshoe distribution
n <- 1000 # number of samples to draw

lambda <- abs(rcauchy(n))
x <- rnorm(n, sd = lambda * tau)
# Or
lambda <- abs(rcauchy(n, scale = tau))
x <- rnorm(n, sd = lambda)

For completeness, this is the R code to reproduce the plot.

set.seed(1234)

# Probability density function of the horseshoe distribution with scale tau.
dhorseshoe <- function(x, tau = 1) {
  theta2 <- x^2 / (2 * tau^2)
  expint::expint_E1(theta2, scale = TRUE) / sqrt(2 * pi^3 * tau^2)
}

tau <- 2 # scale parameter of horseshoe distribution
n <- 1000 # number of samples to draw

lambda <- abs(rcauchy(n))
x <- rnorm(n, sd = lambda * tau)
# Or
lambda <- abs(rcauchy(n, scale = tau))
x <- rnorm(n, sd = lambda)

hist(
  # The horseshoe distribution has heavy tails and is centered at 0.
  # Let's "zoom in" on values between -100 and 100.
  x[abs(x) < 100],
  breaks = 200, probability = TRUE, col = NULL,
  xlim = c(-100, 100), xlab = NULL, ylim = c(0, 0.4), ylab = NULL,
  main = paste("Horseshoe distribution with scale", tau)
)
curve(dhorseshoe(x, tau = tau), n = 1001, add = TRUE, col = "red")