\[ \newcommand{\argmin}{\operatorname{argmin}} \newcommand{\sign}{\operatorname{sign}} \newcommand{\diag}[1]{\operatorname{diag}(#1)} \newcommand{\prox}[2]{\operatorname{prox}_{#1}(#2)} \]

1 Horseshoe Plus

The horseshoe+ prior, as a scale mixture of normals, is given by \[ \begin{align*} (\theta_j | \lambda_j, \tau) &\sim \operatorname{N}(0, \lambda_j^2) \\ \lambda_j &\sim \operatorname{C}^+(0,\eta_j \tau) \\ \eta_j &\sim \operatorname{C}^+(0,1) \end{align*} \] The original horseshoe is \[ \begin{align*} (\theta_j | \lambda_j, \tau) &\sim \operatorname{N}(0, \lambda_j^2) \\ \lambda_j &\sim \operatorname{C}^+(0,\tau) \end{align*} \]

2 Simulations

In the following simulations we generate data according to \[ \begin{align*} y_k &\sim \operatorname{N}\left(\begin{pmatrix}\theta_1 \\ \theta_2\end{pmatrix}, \bf{I}\right) \\ \tau &\sim \operatorname{C}^+(0,1) \end{align*} \]

hist.ci.pct = 0.95
A=10
J=2 
K = 100
theta.true = c(0,0)

We construct a 100 replicates of a signal \(\theta=(0, 0)\).

We generate data with the following:

test.data = list('J'=J, 'K'=K,
                 'y'=t(replicate(K, rnorm(J,theta.true,1))))

2.1 Estimation

First, we load the necessary R packages:

library(rstan)
set_cppo("fast")
library(ggplot2)
library(plyr)
library(reshape2)

The horseshoe+ is implemented in Stan with the following

stan.hsplus.code = "
  data {
    int<lower=0> K; 
    int<lower=0> J; 
    vector[J] y[K]; 
  }
  parameters {
    vector[J] theta_step; 
    vector<lower=0>[J] lambda; 
    vector<lower=0>[J] eta;
    real<lower=0> tau;
  }
  transformed parameters {
    vector[J] theta; 
    theta <- ((theta_step .* lambda) .* eta) * tau;
  }
  model {
    tau ~ cauchy(0, 1);
    eta ~ cauchy(0, 1);
    lambda ~ cauchy(0, 1);
    theta_step ~ normal(0, 1);
    for (k in 1:K) {
      y[k] ~ normal(theta, 1);
    }
  }  
"

and the regular horseshoe:

stan.hs.code = "
  data {
    int<lower=0> K; 
    int<lower=0> J; 
    vector[J] y[K]; 
  }
  parameters {
    vector[J] theta_step; 
    vector<lower=0>[J] lambda; 
    real<lower=0> tau;
  }
  transformed parameters {
    vector[J] theta; 
    theta <- (theta_step .* lambda) * tau;
  }
  model {
    tau ~ cauchy(0, 1);
    lambda ~ cauchy(0, 1);
    theta_step ~ normal(0, 1);
    for (k in 1:K) {
      y[k] ~ normal(theta, 1);
    }
  }  
"

and the \(\operatorname{N}(0,1)\):

stan.norm.code = "
  data {
    int<lower=0> K; 
    int<lower=0> J; 
    vector[J] y[K]; 
  }
  parameters {
    vector[J] theta; 
  }
  model {
    theta ~ normal(0, 1);
    for (k in 1:K) {
      y[k] ~ normal(theta, 1);
    }
  }  
"

It’s necessary to compile the code in Stan (we use Clang):

stan.hsplus.fit = stan_model(model_code=stan.hsplus.code, model_name="hs+ cauchy")
stan.hs.fit = stan_model(model_code=stan.hs.code, model_name="hs cauchy")
stan.norm.fit = stan_model(model_code=stan.norm.code, model_name="normal")
n.iters = 3000
n.chains = 1

Stan is run with 1 chains of 3000 iterations each.

2.2 Horseshoe+ Results

smpls.hsplus.res = sampling(stan.hsplus.fit, 
                            data = test.data, 
                            iter = n.iters,
                            #init = 0,
                            #seed = rng.seed, 
                            chains = n.chains)
theta.smpls.hsplus = extract(smpls.hsplus.res, pars=c("theta"), permuted=TRUE)[[1]]
colnames(theta.smpls.hsplus) = c("theta1", "theta2")

hist.hsplus.ci.pct = 0.95

Next, we produce a density plot of the two coordinates and a histogram of the data within a 95% interval for the product \(\theta_1 \theta_2\) for the prior and posterior distributions:

ggplot(as.data.frame(theta.smpls.hsplus), aes(theta1, theta2)) + geom_density2d()

prior.prod.theta.hsplus = replicate(nrow(theta.smpls.hsplus), 
                                     {
                                     tau = rcauchy(1, 0, 1)
                                     eta = rcauchy(2, 0, 1)
                                     lambda = rcauchy(2, c(0,0), abs(tau * eta))
                                     theta = rnorm(2, lambda, c(1,1)) 
                                     theta[1]*theta[2]
                                     })
post.prod.theta.hsplus = apply(theta.smpls.hsplus, 1, function(x) x[1]*x[2])
prod.dist.hsplus = rbind(data.frame(type="prior", value=prior.prod.theta.hsplus),
                          data.frame(type="posterior", value=post.prod.theta.hsplus))

hist.quants.hsplus = quantile(prod.dist.hsplus$value, probs=c(hist.hsplus.ci.pct, 1.0-hist.hsplus.ci.pct))
hist.data.hsplus = subset(prod.dist.hsplus, min(hist.quants.hsplus) < value & value < max(hist.quants.hsplus))
hist.breaks.hsplus = hist(hist.data.hsplus$value, plot=FALSE, breaks="Scott")$breaks
ggplot(hist.data.hsplus, 
       aes(x=value, group=type)) + geom_histogram(aes(fill=type, y=..density..), breaks=hist.breaks.hsplus) + 
xlab(expression(theta[1]*theta[2]))

The summary statistics for the prior and posterior samples, respectively:

##           Min.        1st Qu.         Median           Mean        3rd Qu. 
## -550581.004800      -2.636474       0.017218    5064.462969       2.581210 
##           Max. 
## 7143408.636000
##             Min.          1st Qu.           Median             Mean 
## -5.040358271e-02 -3.792304437e-04  7.162445000e-08  4.158860680e-05 
##          3rd Qu.             Max. 
##  5.520606967e-04  2.527575922e-02

2.3 Horseshoe Results

smpls.hs.res = sampling(stan.hs.fit, 
                        data = test.data, 
                        iter = n.iters,
                        #init = 0,
                        #seed = rng.seed, 
                        chains = n.chains)
theta.smpls.hs = extract(smpls.hs.res, pars=c("theta"), permuted=TRUE)[[1]]
colnames(theta.smpls.hs) = c("theta1", "theta2")

hist.hs.ci.pct = 0.95
ggplot(as.data.frame(theta.smpls.hs), aes(theta1, theta2)) + geom_density2d()

prior.prod.theta.hs = replicate(nrow(theta.smpls.hs), 
                                     {
                                     tau = rcauchy(1, 0, 1)
                                     lambda = rcauchy(2, c(0,0), abs(tau))
                                     theta = rnorm(2, lambda, c(1,1)) 
                                     theta[1]*theta[2]
                                     })
post.prod.theta.hs = apply(theta.smpls.hs, 1, function(x) x[1]*x[2])
prod.dist.hs = rbind(data.frame(type="prior", value=prior.prod.theta.hs),
                          data.frame(type="posterior", value=post.prod.theta.hs))

hist.quants.hs = quantile(prod.dist.hs$value, probs=c(hist.hs.ci.pct, 1.0-hist.hs.ci.pct))
hist.data.hs = subset(prod.dist.hs, min(hist.quants.hs) < value & value < max(hist.quants.hs))
hist.breaks.hs = hist(hist.data.hs$value, plot=FALSE, breaks="Scott")$breaks
ggplot(hist.data.hs, 
       aes(x=value, group=type)) + geom_histogram(aes(fill=type, y=..density..), breaks=hist.breaks.hs) + 
xlab(expression(theta[1]*theta[2]))