This package provides fast algorithms to calculate the marginal posterior probabilities of data points having a non-zero mean: \[
\pi_n(\theta_i \neq 0 \mid X), \qquad i=1,\ldots,n.
\] Given these, it is easy to calculate the posterior mean, which is also provided: \[
\mathbb{E}[\theta | X] = (\mathbb{E}[\theta_i | X], \ldots, \mathbb{E}[\theta_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 \(\pi_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(n^2)\) 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 \[
\Lambda_n = \text{Beta}(\kappa,\lambda)
\] 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(m n^{3/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 \(\Phi = (\Phi_1,\ldots,\Phi_n)\) and \(\Psi = (\Psi_1, \ldots, \Psi_n)\) that specify the densities of the data points \(X_1,\dots,X_n\) conditional on \(\theta_i\) being \(0\) or \(\theta_i\) being distributed according to \(G\), respectively: \[
\Phi_i = \phi_\sigma(X_i)
\] \[
\Psi_i = \int_{-\infty}^\infty \phi_\sigma(X_i - \theta_i) d G(\theta_i),
\] where \(\phi_\sigma(X_i)\) is the density of a Gaussian with mean \(0\) and standard deviation \(\sigma\). Advanced users may implement their own slab distributions \(G\) by calculating the \(\Psi\) vector themselves. Custom values for \(\Phi\) and \(\Psi\) may also be used to model non-Gaussian noise or to incorporate into \(\Phi\) a shrinkage prior rather than a point-mass at zero.