\[ \newcommand{\argmin}{\operatorname{argmin}} \newcommand{\sign}{\operatorname{sign}} \newcommand{\diag}[1]{\operatorname{diag}(#1)} \newcommand{\prox}[2]{\operatorname{prox}_{#1}(#2)} \]
The Horseshoe+ Prior: Fieller-Creasy
Brandon Willard
University of Chicago Booth17 February, 2015
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 produce 100 replicates of the signal \(\theta=(0, 0)\):
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 = 1Stan 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.95Next, we produce a density plot of the two coordinates and a histogram of the data within a 95% interval for the ratio \(\theta_1/\theta_2\) for the prior and posterior distributions:
ggplot(as.data.frame(theta.smpls.hsplus), aes(theta1, theta2)) + geom_density2d()
prior.ratio.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.ratio.theta.hsplus = apply(theta.smpls.hsplus, 1, function(x) x[1]/x[2])
ratio.dist.hsplus = rbind(data.frame(type="prior", value=prior.ratio.theta.hsplus),
data.frame(type="posterior", value=post.ratio.theta.hsplus))
hist.quants.hsplus = quantile(ratio.dist.hsplus$value, probs=c(hist.hsplus.ci.pct, 1.0-hist.hsplus.ci.pct))
hist.data.hsplus = subset(ratio.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
## -14235.09086000 -1.01917985 -0.00865337 27.65311096
## 3rd Qu. Max.
## 0.85281466 54155.02507000## Min. 1st Qu. Median Mean
## -34535.46763000 -1.40629761 -0.08032056 44.70213712
## 3rd Qu. Max.
## 0.82562027 70537.178760002.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.95ggplot(as.data.frame(theta.smpls.hs), aes(theta1, theta2)) + geom_density2d()
prior.ratio.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.ratio.theta.hs = apply(theta.smpls.hs, 1, function(x) x[1]/x[2])
ratio.dist.hs = rbind(data.frame(type="prior", value=prior.ratio.theta.hs),
data.frame(type="posterior", value=post.ratio.theta.hs))
hist.quants.hs = quantile(ratio.dist.hs$value, probs=c(hist.hs.ci.pct, 1.0-hist.hs.ci.pct))
hist.data.hs = subset(ratio.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]))
The summary statistics for the prior and posterior samples, respectively:
## Min. 1st Qu. Median Mean
## -690.6946057000 -0.8984917765 0.0139350884 0.0933478004
## 3rd Qu. Max.
## 0.9742570663 1767.5718660000## Min. 1st Qu. Median Mean
## -2384.3764460000 -1.3582885170 -0.1262560514 -2.0732908550
## 3rd Qu. Max.
## 0.7163554593 545.58118130003 Normal Results
smpls.norm.res = sampling(stan.norm.fit,
data = test.data,
iter = n.iters,
#init = 0,
#seed = rng.seed,
chains = n.chains)
theta.smpls.norm = extract(smpls.norm.res, pars=c("theta"), permuted=TRUE)[[1]]
colnames(theta.smpls.norm) = c("theta1", "theta2")
hist.norm.ci.pct = 0.95Next, we produce a histogram of the \(\theta\) sum of squares for the prior and posterior distributions:
ggplot(as.data.frame(theta.smpls.norm), aes(theta1, theta2)) + geom_density2d()
prior.ratio.theta.norm = replicate(nrow(theta.smpls.norm),
{
theta = rnorm(2, 0, 1)
theta[1]/theta[2]
})
post.ratio.theta.norm = apply(theta.smpls.norm, 1, function(x) x[1]/x[2])
ratio.dist.norm = rbind(data.frame(type="prior", value=prior.ratio.theta.norm),
data.frame(type="posterior", value=post.ratio.theta.norm))
hist.quants.norm = quantile(ratio.dist.norm$value, probs=c(hist.norm.ci.pct, 1.0-hist.norm.ci.pct))
hist.data.norm = subset(ratio.dist.norm, min(hist.quants.norm) < value & value < max(hist.quants.norm))
hist.breaks.norm = hist(hist.data.norm$value, plot=FALSE, breaks="Scott")$breaks
ggplot(hist.data.norm,
aes(x=value, group=type)) + geom_histogram(aes(fill=type, y=..density..), breaks=hist.breaks.norm) +
xlab(expression(theta[1]/theta[2]))
The summary statistics for the prior and posterior samples, respectively:
## Min. 1st Qu. Median Mean
## -2475.0293020000 -1.0244841590 -0.0096037118 -1.1720561240
## 3rd Qu. Max.
## 0.9124952015 842.2672579000## Min. 1st Qu. Median Mean
## -624.6860726000 -1.6468750710 -0.5716435819 2.0641523620
## 3rd Qu. Max.
## 0.3028431102 2647.0116540000