\[ \newcommand{\argmin}{\operatorname{argmin}} \newcommand{\sign}{\operatorname{sign}} \newcommand{\diag}[1]{\operatorname{diag}(#1)} \newcommand{\prox}[2]{\operatorname{prox}_{#1}(#2)} \]
The Horseshoe+ Prior: Product of Normal Means
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 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 = 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 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-022.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.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]))
The summary statistics for the prior and posterior samples, respectively:
## Min. 1st Qu. Median Mean
## -10868124.950000 -2.068559 -0.022201 -7944.384773
## 3rd Qu. Max.
## 1.802785 300971.074700## Min. 1st Qu. Median Mean
## -4.525801458e-02 -7.048785554e-04 1.787405740e-06 2.149248244e-04
## 3rd Qu. Max.
## 1.079797878e-03 7.607665954e-023 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.prod.theta.norm = replicate(nrow(theta.smpls.norm),
{
theta = rnorm(2, 0, 1)
theta[1]*theta[2]
})
post.prod.theta.norm = apply(theta.smpls.norm, 1, function(x) x[1]*x[2])
prod.dist.norm = rbind(data.frame(type="prior", value=prior.prod.theta.norm),
data.frame(type="posterior", value=post.prod.theta.norm))
hist.quants.norm = quantile(prod.dist.norm$value, probs=c(hist.norm.ci.pct, 1.0-hist.norm.ci.pct))
hist.data.norm = subset(prod.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
## -4.949877701000 -0.400063578400 0.007540719826 -0.024096233750
## 3rd Qu. Max.
## 0.355662903800 5.302094683000## Min. 1st Qu. Median Mean
## -0.0775173415500 -0.0036054260840 0.0005018999935 0.0016326470980
## 3rd Qu. Max.
## 0.0063902109080 0.0640370026200