SequenceSpikeSlab-vignette

library(SequenceSpikeSlab)

Introduction

This package provides fast algorithms for exact Bayesian inference in the sparse normal sequence model. It implements the methods described in (Erven and Szabó 2021). Special care has been taken to make the methods scale to large data sets, and to minimize numerical errors (which arise in all software because floating point numbers are represented with finite precision).

The Sparse Normal Sequence Model

In the sparse normal sequence model we get a sequence of noisy observations X = (X1, …, Xn) such that Xi = θi + ϵi,   i = 1, …, n, where the noise variables ϵi are assumed to be independent with distribution 𝒩(0, σ). We are interested in the underlying means θ = (θ1, …, θn), which are assumed to be sparse. This means that only s of the components in θ are non-zero, with s relatively small compared to n.

The Bayesian Prior

We consider a Bayesian model for the data, with a hierarchical prior on θ of the following form:

  1. First, a sparsity level s ∈ {0, 1, …, n} is chosen according to a prior πn(s).
  2. Then, given s, a subset of non-zero coordinates 𝒮 ⊂ {0, 1, …, n} of size |𝒮| = s is chosen uniformly at random.
  3. Finally, given s and 𝒮, each non-zero mean θi with i ∈ 𝒮 is chosen independently according to a slab prior G, and all other means θi with i ∉ 𝒮 are assumed to be 0.

Given this structure of prior, it thus remains to choose a prior πn on the sparsity level and a slab prior G. We assume that G has a density g. Out of the box, the package supports two choices of slab prior: $$ g_\lambda(\theta_i) = \frac{\lambda}{2} e^{-\lambda |\theta_i|} \qquad \text{(Laplace)} $$ $$ g_\gamma(\theta_i) = \frac{1}{\gamma \pi (1 + (\frac{\theta_i}{\gamma})^2)} \qquad \text{(Cauchy)} $$ Advanced users may also specify their own slab priors.

Theoretically motivated choices for πn (Castillo and Vaart 2012) include πn(s) ∝ e−κslog (3n/s) $$ \pi_n(s) \propto \binom{2n-s}{n}^\kappa $$ for e.g. κ = 0.1 or κ = 1.

Spike-and-Slab Priors

A special case of the general class of priors described above are so-called spike-and-slab priors. Under these priors, there is a mixing hyper-parameter α ∈ [0, 1] and θi|α ∼ (1 − α)δ0 + αG,   i = 1, …, n,

α ∼ Λn. This means that we need to specify a prior Λn on α. Then, conditional on α, the θi are independently distributed. They are 0 with probability 1 − α and they are sampled from the slab distribution G with probability α.

Spike-and-slab priors are a special case of the general hierarchical prior with $$ \pi_n(s) = \binom{n}{s} \int_0^1 \alpha^s (1-\alpha)^{n-s} d \Lambda_n(\alpha). $$

Package Functionality

This package provides fast algorithms to calculate the marginal posterior probabilities of data points having a non-zero mean: πn(θi ≠ 0 ∣ X),   i = 1, …, n. Given these, it is easy to calculate the posterior mean, which is also provided: 𝔼[θ|X] = (𝔼[θi|X], …, 𝔼[θn|X])

Main Methods

There are two main methods: general_sequence_model and fast_spike_slab_beta.

general_sequence_model

This method can handle the general hierarchical prior described above. This means it requires the user to provide a choice for πn as input, and to specify whether the slab prior G should be a Laplace or a Cauchy distribution.

The run time of this method scales as O(n2) in the sample size n. It has been used to handle sample sizes up to 20 000 within half an hour of computation time.

fast_spike_slab_beta

This method is a faster special-purpose method for the spike-and-slab prior with Λn = Beta(κ, λ) This corresponds to the beta-binomial prior $$ \pi_n(s) \propto \frac{\Gamma(\kappa + s) \Gamma(\lambda + n - s)}{s! (n-s)!} $$ in the general hierarchical prior. This method further requires specifying whether the slab prior G should be a Laplace or a Cauchy distribution.

The run time of this method scales as O(mn3/2), where m is a user-supplied parameter that controls the precision of the method. Larger values of m give more precision but slow down the algorithm. The default value of m = 20 has been seen to provide sufficient precision for data sets of different sizes. The method has been used to handle sample sizes up to 100 000 within half and our of computation time.

Advanced Usage

The package is organized as a set of supporting functions in R, which wrap around two C++ functions that implement the main algorithms. The C++ function corresponding to general_sequence_model is accessible directly via SSS_hierarchical_prior and SSS_hierarchical_prior_binomial. The C++ function corresponding to fast_spike_slab_beta can be accessed via SSS_discrete_spike_slab. There are further supporting functions whose name starts with ‘SSS_’.

Implementing Custom Slab Distributions

The C++ functions do not take the data X directly as input. Instead, they require two vectors Φ = (Φ1, …, Φn) and Ψ = (Ψ1, …, Ψn) that specify the densities of the data points X1, …, Xn conditional on θi being 0 or θi being distributed according to G, respectively: Φi = ϕσ(Xi) Ψi = ∫−∞∞ϕσ(Xi − θi)dG(θi), where ϕσ(Xi) is the density of a Gaussian with mean 0 and standard deviation σ. Advanced users may implement their own slab distributions G by calculating the Ψ vector themselves. Custom values for Φ and Ψ may also be used to model non-Gaussian noise or to incorporate into Φ a shrinkage prior rather than a point-mass at zero.

References

Castillo, Ismaël, and Aad van der Vaart. 2012. “Needles and Straw in a Haystack: Posterior Concentration for Possibly Sparse Sequences.” Ann. Statist. 40 (4): 2069–2101. https://doi.org/10.1214/12-AOS1029.
Erven, Tim van, and Botond Szabó. 2021. “Fast Exact Bayesian Inference for Sparse Signals in the Normal Sequence Model.” Bayesian Analysis 16 (3): 933–60. https://doi.org/10.1214/20-BA1227.