HS.MMLE: MMLE for the horseshoe prior for the sparse normal means...
In horseshoe: Implementation of the Horseshoe Prior

Description Usage Arguments Details Value Note References See Also Examples

View source: R/HS_MMLE.R

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).

1	HS.MMLE(y, Sigma2)

`y`	The data, a n1* vector.
`Sigma2`	The variance of the data.

The normal means model is:

y_i=β_i+ε_i, ε_i \sim N(0,σ^2)

And the horseshoe prior:

β_j \sim N(0,σ^2 λ_j^2 τ^2)

λ_j \sim Half-Cauchy(0,1).

This function estimates τ. A plug-in value of σ^2 is used.

The MMLE for the parameter tau of the horseshoe.

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.

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.

The estimated value of τ 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 τ is preferred, see HS.normal.means for the normal means problem, or horseshoe for linear regression.

## 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)