\[ \newcommand{\argmin}{\operatorname{argmin}} \newcommand{\sign}{\operatorname{sign}} \newcommand{\diag}[1]{\operatorname{diag}(#1)} \newcommand{\prox}[2]{\operatorname{prox}_{#1}(#2)} \]
The Horseshoe+ Prior: Normal Vector Mean
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) \end{align*} \] and we consider the following two \(\tau\) priors: \[ \begin{align*} \tau &\sim \operatorname{Unif}[0,1] \\ \tau &\sim \operatorname{C}^+(0,1/J) \end{align*} \]
Qn=10
A=7
J=200
theta.true = rep(0, J)
theta.true[sample.int(J, Qn)] = AWe construct a signal with 10 non-zero components of magnitude 7 and dimension \(J=200\).
We generate the data with the following:
test.data = list('J'=J, 'y'=rnorm(J,theta.true,1))2.1 Estimation
First, we load the necessary R packages:
library(rstan)
library(ggplot2)
theme_set(theme_bw())
library(plyr)
library(dplyr)
library(reshape2)
set_cppo("fast")The horseshoe+ is implemented in Stan with the following
stan.hsplus.unif.code = "
data {
int<lower=1> J;
vector[J] y;
}
parameters {
vector[J] theta_step;
vector[J] lambda_step;
vector[J] eta;
real tau;
}
transformed parameters {
vector[J] theta;
vector[J] lambda;
lambda <- (lambda_step .* eta) * tau;
theta <- ((theta_step .* lambda_step) .* eta) * tau;
}
model {
tau ~ uniform(0, 1);
eta ~ cauchy(0, 1);
lambda_step ~ cauchy(0, 1);
theta_step ~ normal(0, 1);
y ~ normal(theta, 1);
}
"
stan.hsplus.cauchy.code = "
data {
int<lower=1> J;
vector[J] y;
}
parameters {
vector[J] theta_step;
vector[J] lambda_step;
vector[J] eta;
real tau;
}
transformed parameters {
vector[J] theta;
vector[J] lambda;
lambda <- (lambda_step .* eta) * tau;
theta <- ((theta_step .* lambda_step) .* eta) * tau;
}
model {
tau ~ cauchy(0, 1.0/J);
eta ~ cauchy(0, 1);
lambda_step ~ cauchy(0, 1);
theta_step ~ normal(0, 1);
y ~ normal(theta, 1);
}
"and the regular horseshoe:
stan.hs.unif.code = "
data {
int<lower=1> J;
vector[J] y;
}
parameters {
vector[J] theta_step;
vector[J] lambda_step;
real tau;
}
transformed parameters {
vector[J] theta;
vector[J] lambda;
lambda <- lambda_step * tau;
theta <- (theta_step .* lambda_step) * tau;
}
model {
tau ~ uniform(0, 1);
lambda_step ~ cauchy(0, 1);
theta_step ~ normal(0, 1);
y ~ normal(theta, 1);
}
"
stan.hs.cauchy.code = "
data {
int<lower=1> J;
vector[J] y;
}
parameters {
vector[J] theta_step;
vector[J] lambda_step;
real tau;
}
transformed parameters {
vector[J] theta;
vector[J] lambda;
lambda <- lambda_step * tau;
theta <- (theta_step .* lambda_step) * tau;
}
model {
tau ~ cauchy(0, 1.0/J);
lambda_step ~ cauchy(0, 1);
theta_step ~ normal(0, 1);
y ~ normal(theta, 1);
}
"For comparison, we include the DL prior from Section 2.4 of 1:
stan.dl.gamma.code = '
data {
int<lower=1> J;
real<lower=0> A;
vector[J] y;
}
transformed data {
vector<lower=0>[J] Avec;
real<lower=0> JA;
Avec <- rep_vector(A, J);
JA <- J * A;
}
parameters {
vector[J] theta_step;
simplex[J] phi;
vector<lower=0>[J] psi;
real<lower=0> tau;
}
transformed parameters {
vector[J] theta;
for (j in 1:J) {
theta[j] <- ((theta_step[j] .* sqrt(psi[j])) .* phi[j]) * tau;
}
}
model {
tau ~ gamma(JA, 0.5);
psi ~ exponential(0.5);
phi ~ dirichlet(Avec);
theta_step ~ normal(0, 1);
y ~ normal(theta, 1);
}
'It’s necessary to compile the code in Stan (we use Clang):
stan.hsplus.cauchy.fit = stan_model(model_code=stan.hsplus.cauchy.code,
model_name="hs+ cauchy")
stan.hsplus.unif.fit = stan_model(model_code=stan.hsplus.unif.code,
model_name="hs+ unif")stan.hs.cauchy.fit = stan_model(model_code=stan.hs.cauchy.code,
model_name="hs cauchy")
stan.hs.unif.fit = stan_model(model_code=stan.hs.unif.code,
model_name="hs unif")stan.dl.gamma.fit = stan_model(model_code=stan.dl.gamma.code, model_name="dl gamma")n.iters = 1000
n.chains = 3Stan is run with 3 chains of 1000 iterations each.
2.2 Horseshoe+ Results
2.2.1 Cauchy \(\tau\)
smpls.hsplus.res.cauchy = sampling(stan.hsplus.cauchy.fit,
data = test.data,
iter = n.iters,
pars = c("theta", "lambda", "eta", "tau"),
#init = 0,
#seed = rng.seed,
chains = n.chains)
theta.smpls.hsplus.cauchy = extract(smpls.hsplus.res.cauchy, pars=c("theta"), permuted=TRUE)[[1]]
med.theta.hsplus.cauchy = apply(theta.smpls.hsplus.cauchy, 2, median)
hsplus.sample.data.cauchy = melt(extract(smpls.hsplus.res.cauchy, permuted=TRUE))
colnames(hsplus.sample.data.cauchy) = c("iteration", "component", "value", "variable")The following produces plots for each estimated variable:
hsplus.theta.sample.data.cauchy = hsplus.sample.data.cauchy %>%
filter(variable %in% c("theta", "lambda")) %>%
mutate(variable = as.factor(ifelse(variable == "lambda", "kappa", variable)),
value = ifelse(variable == "kappa", 1/(1+abs(value)), value))
hsplus.theta.sample.data.cauchy.2 = hsplus.theta.sample.data.cauchy %>%
group_by(component, variable) %>%
summarise(upper = quantile(value, prob=0.95),
lower = quantile(value, prob=0.05),
middle = mean(value))
hsplus.theta.kappa.cauchy.plot = ggplot(hsplus.theta.sample.data.cauchy.2,
aes(x=component, y=middle, group=component)) + theme_bw() +
geom_pointrange(aes(ymin=lower, ymax=upper), size=0.3) +
facet_grid(variable ~ ., scales="free_y") +
scale_colour_manual(values=c("#D55E00","#0072B2")) +
xlab("") + ylab("")
print(hsplus.theta.kappa.cauchy.plot)
ggplot(hsplus.sample.data.cauchy %>% filter(variable == "eta") %>%
mutate(value = abs(value)),
aes(x=component, y=value, group=component)) +
geom_boxplot(outlier.size=0.5, outlier.colour="grey") +
scale_y_log10() +
geom_hline(aes(yintercept=0), size=0.5, linetype="dashed") +
geom_hline(aes(yintercept=1), size=0.5, linetype="dashed") +
xlab(expression(eta[i])) + ylab("")
The median \(\theta\) value for the horseshoe+ fit has an SSE of 23.9152969919731.
hsplus.hist.plot.subset.cauchy = hsplus.sample.data.cauchy %>%
filter(variable == "tau") %>%
mutate(value = abs(value))
hsplus.hist.breaks.cauchy = hist(hsplus.hist.plot.subset.cauchy$value,
plot=FALSE, breaks="FD")$breaks
ggplot(hsplus.hist.plot.subset.cauchy, aes(x=value)) +
geom_histogram(breaks = hsplus.hist.breaks.cauchy) + xlab(expression(tau)) 
2.2.2 Uniform \(\tau\)
smpls.hsplus.res.unif = sampling(stan.hsplus.unif.fit,
data = test.data,
iter = n.iters,
pars = c("theta", "lambda", "eta", "tau"),
#init = 0,
#seed = rng.seed,
chains = n.chains)
theta.smpls.hsplus.unif = extract(smpls.hsplus.res.unif, pars=c("theta"), permuted=TRUE)[[1]]
med.theta.hsplus.unif = apply(theta.smpls.hsplus.unif, 2, median)
hsplus.sample.data.unif = melt(extract(smpls.hsplus.res.unif, permuted=TRUE))
colnames(hsplus.sample.data.unif) = c("iteration", "component", "value", "variable")The following produces plots for each estimated variable:
hsplus.theta.sample.data.unif = hsplus.sample.data.unif %>%
filter(variable %in% c("theta", "lambda")) %>%
mutate(variable = as.factor(ifelse(variable == "lambda", "kappa", variable)),
value = ifelse(variable == "kappa", 1/(1+abs(value)), value))
hsplus.theta.sample.data.unif.2 = hsplus.theta.sample.data.unif %>%
group_by(component, variable) %>%
summarise(upper = quantile(value, prob=0.95),
lower = quantile(value, prob=0.05),
middle = mean(value))
hsplus.theta.kappa.plot = ggplot(hsplus.theta.sample.data.unif.2,
aes(x=component, y=middle, group=component)) + theme_bw() +
geom_pointrange(aes(ymin=lower, ymax=upper), size=0.3) +
facet_grid(variable ~ ., scales="free_y") +
scale_colour_manual(values=c("#D55E00","#0072B2")) +
xlab("") + ylab("")
print(hsplus.theta.kappa.plot)
ggplot(hsplus.sample.data.unif %>%
filter(variable == "eta") %>%
mutate(value = abs(value)),
aes(x=component, y=value, group=component)) +
geom_boxplot(outlier.size=0.5, outlier.colour="grey") + scale_y_log10() +
geom_hline(aes(yintercept=0), size=0.5, linetype="dashed") +
geom_hline(aes(yintercept=1), size=0.5, linetype="dashed") +
xlab(expression(eta[i])) + ylab("")
The median \(\theta\) value for the horseshoe+ estimates has an SSE of 23.9264802839268.
hsplus.hist.plot.subset.unif = hsplus.sample.data.unif %>%
filter(variable == "tau") %>%
mutate(value = abs(value))
hsplus.hist.breaks.unif = hist(hsplus.hist.plot.subset.unif$value,
plot=FALSE, breaks="FD")$breaks
ggplot(hsplus.hist.plot.subset.unif, aes(x=value)) +
geom_histogram(breaks = hsplus.hist.breaks.unif) + xlab(expression(tau)) 
2.3 Horseshoe Results
2.3.1 Cauchy \(\tau\)
smpls.hs.res.cauchy = sampling(stan.hs.cauchy.fit,
data = test.data,
iter = n.iters,
pars = c("theta", "lambda", "tau"),
#init = 0,
#seed = rng.seed,
chains = n.chains)
theta.smpls.hs.cauchy = extract(smpls.hs.res.cauchy, pars=c("theta"), permuted=TRUE)[[1]]
med.theta.hs.cauchy = apply(theta.smpls.hs.cauchy, 2, median)
hs.sample.data.cauchy = melt(extract(smpls.hs.res.cauchy, permuted=TRUE))
colnames(hs.sample.data.cauchy) = c("iteration", "component", "value", "variable")The following produces plots for each estimated variable:
hs.theta.sample.data.cauchy = hs.sample.data.cauchy %>%
filter(variable %in% c("theta", "lambda")) %>%
mutate(variable = as.factor(ifelse(variable == "lambda", "kappa", variable)),
value = ifelse(variable == "kappa", 1/(1+abs(value)), value))
hs.theta.sample.data.cauchy.2 = hs.theta.sample.data.cauchy %>%
group_by(component, variable) %>%
summarise(upper = quantile(value, prob=0.95),
lower = quantile(value, prob=0.05),
middle = mean(value))
hs.theta.kappa.cauchy.plot = ggplot(hs.theta.sample.data.cauchy.2,
aes(x=component, y=middle, group=component)) + theme_bw() +
geom_pointrange(aes(ymin=lower, ymax=upper), size=0.3) +
facet_grid(variable ~ ., scales="free_y") +
scale_colour_manual(values=c("#D55E00","#0072B2")) +
xlab("") + ylab("")
print(hs.theta.kappa.cauchy.plot)
The median \(\theta\) value for the horseshoe estimates has an SSE of 23.0325396977138.
hs.hist.plot.subset.cauchy = hs.sample.data.cauchy %>%
filter(variable == "tau") %>%
mutate(value = abs(value))
hs.hist.breaks.cauchy = hist(hs.hist.plot.subset.cauchy$value,
plot=FALSE, breaks="FD")$breaks
ggplot(hs.hist.plot.subset.cauchy, aes(x=value)) +
geom_histogram(breaks = hs.hist.breaks.cauchy) + xlab(expression(tau)) 
2.3.2 Uniform \(\tau\)
smpls.hs.res.unif = sampling(stan.hs.unif.fit,
data = test.data,
iter = n.iters,
pars = c("theta", "lambda", "tau"),
#init = 0,
#seed = rng.seed,
chains = n.chains)
theta.smpls.hs.unif = extract(smpls.hs.res.unif, pars=c("theta"), permuted=TRUE)[[1]]
med.theta.hs.unif = apply(theta.smpls.hs.unif, 2, median)
hs.sample.data.unif = melt(extract(smpls.hs.res.unif, permuted=TRUE))
colnames(hs.sample.data.unif) = c("iteration", "component", "value", "variable")hs.theta.sample.data.unif = hs.sample.data.unif %>%
filter(variable %in% c("theta", "lambda")) %>%
mutate(variable = as.factor(ifelse(variable == "lambda", "kappa", variable)),
value = ifelse(variable == "kappa", 1/(1+abs(value)), value))
hs.theta.sample.data.unif.2 = hs.theta.sample.data.unif %>%
group_by(component, variable) %>%
summarise(upper = quantile(value, prob=0.95),
lower = quantile(value, prob=0.05),
middle = mean(value))
hs.theta.kappa.unif.plot = ggplot(hs.theta.sample.data.unif.2,
aes(x=component, y=middle, group=component)) + theme_bw() +
geom_pointrange(aes(ymin=lower, ymax=upper), size=0.3) +
facet_grid(variable ~ ., scales="free_y") +
scale_colour_manual(values=c("#D55E00","#0072B2")) +
xlab("") + ylab("")
print(hs.theta.kappa.unif.plot)
The median \(\theta\) value for the horseshoe estimates has an SSE of 23.8932365735653.
hs.hist.plot.subset.unif = hs.sample.data.unif %>%
filter(variable == "tau") %>%
mutate(value = abs(value))
hs.hist.breaks.unif = hist(hs.hist.plot.subset.unif$value,
plot=FALSE, breaks="FD")$breaks
ggplot(hs.hist.plot.subset.unif, aes(x=value)) +
geom_histogram(breaks = hs.hist.breaks.unif) + xlab(expression(tau))
2.4 DL Results
#dl.n.iters = 10000
dl.test.data = test.data
dl.test.data$A = 1/J
smpls.dl.res.gamma = sampling(stan.dl.gamma.fit,
data = dl.test.data,
iter = n.iters, #dl.n.iters,
pars = c("theta", "psi", "phi", "tau"),
#init = 0,
#seed = rng.seed,
chains = n.chains)
theta.smpls.dl.gamma = extract(smpls.dl.res.gamma, pars=c("theta"), permuted=TRUE)[[1]]
med.theta.dl.gamma = apply(theta.smpls.dl.gamma, 2, median)
dl.sample.data.gamma = melt(extract(smpls.dl.res.gamma, permuted=TRUE))
colnames(dl.sample.data.gamma) = c("iteration", "component", "value", "variable")dl.theta.sample.data.gamma = dl.sample.data.gamma %>%
filter(variable %in% c("theta", "psi", "phi"))
ggplot(dl.theta.sample.data.gamma,
aes(x=component, y=value, group=component)) +
geom_boxplot(outlier.size=0.5, outlier.colour="grey") +
facet_grid(variable ~ ., scales="free_y") +
geom_point(data=data.frame(value=theta.true,
component=seq_along(theta.true),
variable="theta",
type="true"),
aes(colour=type),
size=1) +
geom_point(data=data.frame(value=drop(test.data$y),
variable="theta",
component=seq_along(drop(test.data$y)),
type="obs"),
aes(colour=type),
size=1) +
scale_colour_manual(values=c("#D55E00","#0072B2")) +
xlab("i") + ylab("")
The median \(\theta\) value for the DL prior estimates has an SSE of 21.9216130035049.
dl.hist.plot.subset.gamma = subset(dl.sample.data.gamma, variable %in% "tau")
dl.hist.breaks.gamma = hist(dl.hist.plot.subset.gamma$value,
plot=FALSE, breaks="FD")$breaks
ggplot(dl.hist.plot.subset.gamma, aes(x=value)) +
geom_histogram(breaks = dl.hist.breaks.gamma) + xlab(expression(tau))
3 Replicated Simulations
We now proceed to replicate simulations over an array of similar problems:
Ns = c(100, 200)
q.dns = c(.01,.05,.1,.2,.3)
As = c(7,8)
rep.reps = 100
rep.iters = 5000
dl.rep.iters = 6000
rep.grid = expand.grid(N=Ns, q.dn=q.dns, A=As)
rep.ex.data = alply(rep.grid, 1,
function(d) {
qq = ceiling(d$N * d$q.dn)
theta.true = rep(0, d$N)
theta.true[sample.int(d$N, qq)] = d$A
J = d$N
return(replicate(rep.reps,
list('A'=d$A,'n'=d$N,'q'=qq,
'theta.true'=theta.true,
'data'=list('J'=J, 'A'=1/J,
'y'=rnorm(J, theta.true, 1))
)
, simplify=FALSE))
})
compute.replications = function(model, data,
rep.iters = 10000, n.chains = 1) {
rep.results = ldply(data,
function(d) {
mcmc.reps.res = ldply(d,
function(dd) {
smpls.res = sampling(model,
data = dd$data,
iter = rep.iters,
pars = c("theta"),
#init = 0,
chains = n.chains)
n.eff = summary(smpls.res)$summary[,"n_eff"]
smpls = extract(smpls.res, pars=c("theta"),
permuted=TRUE)[[1]]
med.theta = apply(smpls, 2, median)
med.err.vals = dd$theta.true - med.theta
SSE.median = crossprod(med.err.vals)
mean.theta = apply(smpls, 2, median)
mean.err.vals = dd$theta.true - mean.theta
SSE.mean = crossprod(mean.err.vals)
return(data.frame("SSE.median"=SSE.median,
"SSE.mean"=SSE.mean,
"N.eff.mean"=mean(n.eff),
"N.eff.min"=min(n.eff)))
}, .parallel=TRUE)
SSE.mean.ci = quantile(mcmc.reps.res$SSE.mean,
probs=c(0.05, 0.5, 0.95))
SSE.median.ci = quantile(mcmc.reps.res$SSE.median,
probs=c(0.05, 0.5, 0.95))
res = data.frame("A"=d[[1]]$A, "n"=d[[1]]$n, "q"=d[[1]]$q,
"N.eff.mean"=mean(mcmc.reps.res$N.eff.mean),
"N.eff.min.mean"=mean(mcmc.reps.res$N.eff.mean),
"SSE.mean.mean"=mean(mcmc.reps.res$SSE.mean),
"SSE.mean.ci.05"=SSE.mean.ci[1],
"SSE.mean.ci.50"=SSE.mean.ci[2],
"SSE.mean.ci.95"=SSE.mean.ci[3],
"SSE.median.mean"=mean(mcmc.reps.res$SSE.median),
"SSE.median.ci.05"=SSE.median.ci[1],
"SSE.median.ci.50"=SSE.median.ci[2],
"SSE.median.ci.95"=SSE.median.ci[3]
)
return(res)
})
return(rep.results)
}resulting in the following subproblems, each replicated 100 times:
rep.grid## N q.dn A
## 1 100 0.01 7
## 2 200 0.01 7
## 3 100 0.05 7
## 4 200 0.05 7
## 5 100 0.10 7
## 6 200 0.10 7
## 7 100 0.20 7
## 8 200 0.20 7
## 9 100 0.30 7
## 10 200 0.30 7
## 11 100 0.01 8
## 12 200 0.01 8
## 13 100 0.05 8
## 14 200 0.05 8
## 15 100 0.10 8
## 16 200 0.10 8
## 17 100 0.20 8
## 18 200 0.20 8
## 19 100 0.30 8
## 20 200 0.30 8We estimate each problem with Stan in parallel and some setup is required:
library(doMC)
registerDoMC(cores=detectCores())3.1 Priors
3.1.1 Uniform \(\tau\)
reps.hsplus.unif.time = system.time({
theta.ss.hsplus.unif = compute.replications(stan.hsplus.unif.fit, rep.ex.data, rep.iters=rep.iters)
})
save(theta.ss.hsplus.unif, file=paste0("theta-ss-hsplus-unif-", rep.iters, ".RData"))reps.hs.unif.time = system.time({
theta.ss.hs.unif = compute.replications(stan.hs.unif.fit, rep.ex.data, rep.iters=rep.iters)
})
save(theta.ss.hs.unif, file=paste0("theta-ss-hs-unif-", rep.iters, ".RData"))rep.hs.unif.mse.data = rbind(cbind(prior="hs+ unif", theta.ss.hsplus.unif),
cbind(prior="hs unif", theta.ss.hs.unif))
rep.hs.unif.mse.data = transform(rep.hs.unif.mse.data,
N=as.factor(N), q.pct=as.factor(q.dn),
A=as.factor(A), prior=as.factor(prior)) 3.1.2 Cauchy \(\tau\)
reps.hsplus.cauchy.time = system.time({
theta.ss.hsplus.cauchy = compute.replications(stan.hsplus.cauchy.fit, rep.ex.data, rep.iters=rep.iters)
})
save(theta.ss.hsplus.cauchy, file=paste0("theta-ss-hsplus-cauchy-", rep.iters, ".RData"))reps.hs.cauchy.time = system.time({
theta.ss.hs.cauchy = compute.replications(stan.hs.cauchy.fit, rep.ex.data, rep.iters=rep.iters)
})
save(theta.ss.hs.cauchy, file=paste0("theta-ss-hs-cauchy-", rep.iters, ".RData"))rep.hs.cauchy.mse.data = rbind(cbind(prior="hs+ cauchy", theta.ss.hsplus.cauchy),
cbind(prior="hs cauchy", theta.ss.hs.cauchy))
rep.hs.cauchy.mse.data = transform(rep.hs.cauchy.mse.data,
N=as.factor(N), q.pct=as.factor(q.dn),
A=as.factor(A), prior=as.factor(prior)) 3.2 Results
Finally, we combine the results for both models to produce a table (in order to get the xtabs table to render properly, the development version of pander on Github needed to be installed):
library(pander)
rep.mse.data = rbind(rep.hs.cauchy.mse.data,
rep.hs.unif.mse.data)
library(tables)
rep.mse.table = tabular((N*q.pct*A*(Heading()*SSE.median.mean))*Heading()*identity ~ prior,
data=rep.mse.data)
Hmisc::html(rep.mse.table, file=stdout())| prior | ||||||
|---|---|---|---|---|---|---|
| N | q.pct | A | hs cauchy | hs+ cauchy | hs unif | hs+ unif |
| 100 | 0.01 | 7 | 1.314 | 1.368 | 1.609 | 2.011 |
| 8 | 1.098 | 1.234 | 1.472 | 2.647 | ||
| 0.05 | 7 | 7.549 | 8.462 | 8.641 | 9.264 | |
| 8 | 7.114 | 7.651 | 7.956 | 7.900 | ||
| 0.1 | 7 | 16.715 | 16.123 | 18.376 | 17.168 | |
| 8 | 14.987 | 15.887 | 16.555 | 15.736 | ||
| 0.2 | 7 | 33.987 | 29.543 | 32.899 | 30.138 | |
| 8 | 35.915 | 33.156 | 35.082 | 32.673 | ||
| 0.3 | 7 | 53.803 | 46.433 | 47.494 | 46.454 | |
| 8 | 52.134 | 43.402 | 44.827 | 41.851 | ||
| 200 | 0.01 | 7 | 2.915 | 3.703 | 3.517 | 4.876 |
| 8 | 2.478 | 3.577 | 2.967 | 3.891 | ||
| 0.05 | 7 | 15.518 | 16.202 | 16.415 | 17.118 | |
| 8 | 14.548 | 14.184 | 15.576 | 16.601 | ||
| 0.1 | 7 | 32.430 | 28.922 | 33.947 | 33.575 | |
| 8 | 33.731 | 31.449 | 35.618 | 33.206 | ||
| 0.2 | 7 | 73.370 | 63.879 | 71.546 | 64.799 | |
| 8 | 67.801 | 57.716 | 65.871 | 61.394 | ||
| 0.3 | 7 | 105.259 | 89.753 | 94.581 | 85.955 | |
| 8 | 103.900 | 88.112 | 105.055 | 84.512 | ||
“Dirichlet-Laplace priors for optimal shrinkage”, Anirban Bhattacharya, Debdeep Pati, Natesh S. Pillai, David B. Dunson↩