HS.MMLE {horseshoe}

MMLE for the horseshoe prior for the sparse normal means problem.

Description

Compute the marginal maximum likelihood estimator (MMLE) of tau for the horseshoe for the normal means problem (i.e. linear regression with the design matrix equal to the identity matrix). The MMLE is explained and studied in Van der Pas et al. (2016).

Usage

HS.MMLE(y, Sigma2)

Arguments

y

The data, a n*1 vector.

Sigma2

The variance of the data.

Details

The normal means model is:

y_i=\beta_i+\epsilon_i, \epsilon_i \sim N(0,\sigma^2)

And the horseshoe prior:

\beta_j \sim N(0,\sigma^2 \lambda_j^2 \tau^2)

\lambda_j \sim Half-Cauchy(0,1).

This function estimates \tau. A plug-in value of \sigma^2 is used.

Value

The MMLE for the parameter tau of the horseshoe.

Note

Requires a minimum of 2 observations. May return an error for vectors of length larger than 400 if the truth is very sparse. In that case, try HS.normal.means.

References

van der Pas, S.L., Szabo, B., and van der Vaart, A. (2017), Uncertainty quantification for the horseshoe (with discussion). Bayesian Analysis 12(4), 1221-1274.

van der Pas, S.L., Szabo, B., and van der Vaart A. (2017), Adaptive posterior contraction rates for the horseshoe. Electronic Journal of Statistics 10(1), 3196-3225.

See Also

The estimated value of \tau can be plugged into HS.post.mean to obtain the posterior mean, and into HS.post.var to obtain the posterior variance. These functions are all for empirical Bayes; if a full Bayes version with a hyperprior on \tau is preferred, see HS.normal.means for the normal means problem, or horseshoe for linear regression.

Examples

## Not run: #Example with 5 signals, rest is noise
truth <- c(rep(0, 95), rep(8, 5))
y <-  truth + rnorm(100)
(tau.hat <- HS.MMLE(y, 1)) #returns estimate of tau
plot(y, HS.post.mean(y, tau.hat, 1)) #plot estimates against the data

## End(Not run)
## Not run: #Example where the data variance is estimated first
truth <- c(rep(0, 950), rep(8, 50))
y <-  truth + rnorm(100, mean = 0, sd = sqrt(2))
sigma2.hat <- var(y)
(tau.hat <- HS.MMLE(y, sigma2.hat)) #returns estimate of tau
plot(y, HS.post.mean(y, tau.hat, sigma2.hat)) #plot estimates against the data

## End(Not run)