Metadata-Version: 2.4
Name: sparsevb
Version: 0.1.3
Summary: Spike-and-Slab Variational Bayes for Linear and Logistic Regression
Author-email: Gabriel Clara <gabriel.j.clara@gmail.com>
License-Expression: GPL-3.0-or-later
Project-URL: Repository, https://gitlab.com/gclara/varpack
Keywords: variational-bayes,spike-and-slab,sparse-regression,bayesian,variable-selection
Classifier: Development Status :: 3 - Alpha
Classifier: Intended Audience :: Science/Research
Classifier: Topic :: Scientific/Engineering :: Mathematics
Classifier: Programming Language :: Python :: 3
Classifier: Programming Language :: Python :: 3.8
Classifier: Programming Language :: Python :: 3.9
Classifier: Programming Language :: Python :: 3.10
Classifier: Programming Language :: Python :: 3.11
Classifier: Programming Language :: Python :: 3.12
Classifier: Programming Language :: Python :: 3.13
Classifier: Programming Language :: Cython
Classifier: Programming Language :: C++
Classifier: Operating System :: POSIX :: Linux
Classifier: Operating System :: MacOS
Requires-Python: >=3.8
Description-Content-Type: text/markdown
Requires-Dist: numpy>=1.20.0
Requires-Dist: scipy>=1.6.0
Requires-Dist: scikit-learn>=1.0.0

# sparsevb

A Python package for **Spike-and-Slab Variational Bayes** for sparse variable selection in linear and logistic regression.

## Description

This package implements the mean-field variational Bayes algorithm of [Ray and Szabo (2022)](https://doi.org/10.1080/01621459.2020.1847121) for scalable Bayesian variable selection. Given a potentially high-dimensional design matrix *X* and response *y*, the algorithm approximates the posterior distribution under a spike-and-slab prior using coordinate-ascent variational inference (CAVI).

The spike-and-slab prior on each coefficient takes the form:

```
theta_i ~ gamma_i * Slab(lambda) + (1 - gamma_i) * delta_0
```

where `gamma_i ~ Bernoulli(r)` is the inclusion indicator and the slab is a Laplace (default) or Gaussian distribution. The variational approximation returns posterior inclusion probabilities `gamma_i` for each variable, along with posterior mean (`mu`) and standard deviation (`sigma`) parameters.

It is a direct port of the [R `sparsevb` package](https://cran.r-project.org/package=sparsevb), using the same underlying C++ implementation with a Pythonic API.

## Installation

```bash
pip install sparsevb
```

**System requirements**: LAPACK, BLAS, and a C++14 compiler (usually available by default on Linux and macOS).

To install from source:

```bash
pip install .
```

## Quick start

### Linear regression

```python
import numpy as np
from sparsevb import svb_fit_linear

# Synthetic data: 200 samples, 20 features, 3 active
np.random.seed(42)
n, p = 200, 20
X = np.random.randn(n, p)
beta_true = np.zeros(p)
beta_true[:3] = [1.5, -2.0, 3.0]
y = X @ beta_true + 0.5 * np.random.randn(n)

# Fit sparse model
res = svb_fit_linear(X, y)

# Point estimate of coefficients
print("Estimated coefficients:", res['mu'] * res['gamma'])

# Variable selection: gamma_i > 0.5 indicates inclusion
selected = res['gamma'] > 0.5
print("Selected variables:", np.where(selected)[0])
```

### Logistic regression

```python
from sparsevb import svb_fit_logit

# Binary response
logits = X @ beta_true
probs = 1 / (1 + np.exp(-logits))
y_binary = (np.random.rand(n) < probs).astype(float)

res = svb_fit_logit(X, y_binary)
print("Estimated coefficients:", res['mu'] * res['gamma'])
```

### With intercept

```python
y_with_intercept = X @ beta_true + 5.0 + 0.5 * np.random.randn(n)
res = svb_fit_linear(X, y_with_intercept, intercept=True)
print("Intercept:", res['intercept'])
print("Coefficients:", res['mu'] * res['gamma'])
```

## API Reference

### `svb_fit_linear(X, y, **kwargs)`

Fit a sparse linear regression model `Y = X @ theta + noise`.

**Parameters:**

| Parameter | Type | Default | Description |
|-----------|------|---------|-------------|
| `X` | array (n, p) | *required* | Design matrix |
| `y` | array (n,) | *required* | Response vector |
| `max_iter` | int | 1000 | Maximum CAVI iterations |
| `tol` | float | 1e-5 | Convergence tolerance on the max change in binary entropy of gamma |
| `slab` | str | `"laplace"` | Slab distribution: `"laplace"` or `"gaussian"`. Laplace is recommended |
| `intercept` | bool | False | Include an intercept term |
| `mu` | array (p,) | None | Initial mean parameters. If None, initialized via Ridge regression |
| `sigma` | array (p,) | None | Initial std parameters. If None, set to ones |
| `gamma` | array (p,) | None | Initial inclusion probabilities. If None, initialized via Lasso |
| `noise_sd` | float | None | Known noise std. If None, estimated from data |

**Returns** a dict with keys:

| Key | Type | Description |
|-----|------|-------------|
| `mu` | ndarray (p,) | Posterior mean of the slab component |
| `sigma` | ndarray (p,) | Posterior std of the slab component |
| `gamma` | ndarray (p,) | Posterior inclusion probabilities in [0, 1] |
| `noise_sd` | float | Noise std used (estimated or provided) |
| `intercept` | float or None | Estimated intercept value |

### `svb_fit_logit(X, y, **kwargs)`

Fit a sparse logistic regression model `P(y=1 | X) = logistic(X @ theta)`.

Same parameters as `svb_fit_linear` except:
- `y` must be binary (0 or 1)
- `noise_sd` is not available (no noise parameter in logistic regression)

Returns the same dict, with `noise_sd` always set to None.

## Interpreting the output

- **Point estimate**: the posterior mean of theta under the variational approximation is `mu * gamma`.
- **Variable selection**: a variable *i* is typically selected if `gamma[i] > 0.5`.
- **Posterior uncertainty**: `sigma[i]` gives the posterior standard deviation of the slab component for coefficient *i*. Note that variational Bayes tends to underestimate posterior variance.

## Algorithm details

The algorithm uses **coordinate-ascent variational inference** (CAVI) to minimize the KL divergence between the mean-field variational family and the true posterior. Each iteration updates `mu_i`, `sigma_i`, and `gamma_i` for every coordinate (see Algorithm 1 in Ray & Szabo, 2022).

Key implementation choices:

- **Laplace slab prior**: unlike many VB approaches that use Gaussian slabs, this implementation defaults to Laplace slabs. Laplace priors avoid the excessive shrinkage of non-zero coefficients that Gaussian slabs can cause, leading to better signal recovery.
- **Prioritized updating order**: the CAVI algorithm is sensitive to the order in which coordinates are updated. This implementation sorts coordinates by decreasing `|mu_init|` (from a preliminary Ridge estimate), updating the most important coefficients first. This avoids poor local minima that can arise with lexicographic or random orderings.
- **Noise estimation**: for linear regression with unknown noise variance, the data is rescaled by an estimate of the noise standard deviation (empirical Bayes approach). When `n > p`, Ridge residuals are used; otherwise, the MAD of Lasso residuals is used.
- **Convergence**: the algorithm stops when the maximum change in binary entropy of the inclusion probabilities falls below `tol`.

## Dependencies

- Python >= 3.8
- NumPy >= 1.20.0
- SciPy >= 1.6.0
- scikit-learn >= 1.0.0

Build dependencies: Cython >= 0.29.0, a C++14 compiler, LAPACK, BLAS.

## License

GPL-3.0-or-later

## References

- Ray, K. and Szabo, B. (2022). **Variational Bayes for high-dimensional linear regression with sparse priors.** *Journal of the American Statistical Association*, 117(539), 1270-1281. [DOI:10.1080/01621459.2020.1847121](https://doi.org/10.1080/01621459.2020.1847121)
