\[ \newcommand{\argmin}{\operatorname{argmin}} \newcommand{\sign}{\operatorname{sign}} \newcommand{\diag}[1]{\operatorname{diag}(#1)} \newcommand{\prox}[2]{\operatorname{prox}_{#1}(#2)} \]
The Horseshoe+ Prior: Efron Maximum Problem
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(\theta, \bf{I}\right) \\ \tau &\sim \operatorname{C}^+(0,1) \end{align*} \]
Qn=1
A=10
J=100
theta.true = c(rep(A,Qn),rep(0,J-Qn))We construct a signal with 1 non-zero components of magnitude 10 and dimension \(J=100\).
We generate the data with the following:
test.data = list('J'=J, 'K'=1,
'y'=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> J;
vector[J] y;
}
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);
y ~ normal(theta, 1);
}
"and the regular horseshoe:
stan.hs.code = "
data {
int<lower=0> J;
vector[J] y;
}
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);
y ~ normal(theta, 1);
}
"and the \(\operatorname{N}(0,300)\):
stan.norm.code = "
data {
int<lower=0> J;
vector[J] y;
}
parameters {
vector[J] theta;
}
model {
theta ~ normal(0, sqrt(300));
y ~ 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 = 1500
n.chains = 1Stan is run with 1 chains of 1500 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]]
hsplus.sample.data = melt(extract(smpls.hsplus.res, permuted=TRUE))
colnames(hsplus.sample.data) = c("iteration", "component", "value", "variable")
hist.hsplus.ci.pct = 0.60Next, we produce a histogram of the data within a 60% interval of the \(\theta\) sum of squares for the prior and posterior distributions:
post.max.theta.hsplus = apply(theta.smpls.hsplus, 1, function(x) max(x))
prior.max.theta.hsplus = replicate(nrow(theta.smpls.hsplus),
{
tau = rcauchy(1, 0, 1)
eta = rcauchy(J, 0, 1)
lambda = rcauchy(J, 0, abs(tau * eta))
theta = rnorm(J, lambda, 1)
return(max(theta))
})
max.dist.hsplus = rbind(data.frame(type="prior", value=prior.max.theta.hsplus),
data.frame(type="posterior", value=post.max.theta.hsplus))
hist.quants.hsplus = quantile(max.dist.hsplus$value, probs=c(hist.hsplus.ci.pct, 1.0-hist.hsplus.ci.pct))
hist.data.hsplus = subset(max.dist.hsplus, 0 < 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("max")
The summary statistics for the prior and posterior max samples, respectively:
## Min. 1st Qu. Median Mean
## 1.7492655 53.1268591 200.8580190 6541.0574080
## 3rd Qu. Max.
## 755.8247259 1549310.4860000## sd=75279.342402## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 7.93647013 10.44591069 11.13987144 11.11395183 11.76093402 14.45372138## sd=1.0072232.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]]
hs.sample.data = melt(extract(smpls.hs.res, permuted=TRUE))
colnames(hs.sample.data) = c("iteration", "component", "value", "variable")
hist.hs.ci.pct = 0.70Next, we produce a histogram of the data within a 70% interval of the \(\theta\) sum of squares for the prior and posterior distributions:
post.max.theta.hs = apply(theta.smpls.hs, 1, function(x) max(x))
prior.max.theta.hs = replicate(nrow(theta.smpls.hs),
{
tau = rcauchy(1, 0, 1)
lambda = rcauchy(J, 0, abs(tau))
theta = rnorm(J, lambda, 1)
return(max(theta))
})
max.dist.hs = rbind(data.frame(type="prior", value=prior.max.theta.hs),
data.frame(type="posterior", value=post.max.theta.hs))
hist.quants.hs = quantile(max.dist.hs$value, probs=c(hist.hs.ci.pct, 1.0-hist.hs.ci.pct))
hist.data.hs = subset(max.dist.hs, 0 < 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("max")
The summary statistics for the prior and posterior max samples, respectively:
## Min. 1st Qu. Median Mean 3rd Qu.
## 2.074140 15.331369 51.663985 7218.170104 183.214979
## Max.
## 4665468.160000## sd=170506.003884## Min. 1st Qu. Median Mean 3rd Qu.
## 7.565786932 10.459012860 11.113216320 11.084773380 11.734640740
## Max.
## 14.819003870## sd=1.0479393 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]]
norm.sample.data = melt(extract(smpls.norm.res, permuted=TRUE))
colnames(norm.sample.data) = c("iteration", "component", "value", "variable")Next, we produce a histogram of the \(\theta\) sum of squares for the prior and posterior distributions:
post.max.theta.norm = apply(theta.smpls.norm, 1, function(x) max(x))
prior.max.theta.norm = replicate(nrow(theta.smpls.norm),
{
theta = rnorm(J, 0, 1)
return(max(theta))
})
max.dist.norm = rbind(data.frame(type="prior", value=prior.max.theta.norm),
data.frame(type="posterior", value=post.max.theta.norm))
hist.data.norm = max.dist.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("max")
The summary statistics for the prior and posterior max samples, respectively:
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 1.275850788 2.245103688 2.478996202 2.530693836 2.771210479 4.172248284## sd=0.412846## Min. 1st Qu. Median Mean 3rd Qu.
## 8.306781272 10.560446410 11.265102400 11.285748420 11.995709810
## Max.
## 14.605943400## sd=1.018226