HorseshoePrior#

class cultivars.bayes.priors.HorseshoePrior(*, deterministic_scale=10.0)[source]#

The horseshoe: aggressive shrinkage that lets large signals through.

Each standardized lag coefficient gets a half-Cauchy local scale and the whole lag block shares a half-Cauchy global one (Carvalho, Polson & Scott 2010),

\[\beta_{ij} \mid \lambda_{ij}, \tau \sim \mathcal{N}(0, \tau^2 \lambda_{ij}^2), \qquad \lambda_{ij} \sim \mathcal{C}^+(0, 1), \qquad \tau \sim \mathcal{C}^+(0, 1),\]

where \(\beta_{ij}\) is the coefficient divided by the Minnesota unit ratio \(s_i / s_j\), so one latent scale means the same thing whatever the variables’ units. The Cauchy tails are the substance. The shrinkage weight \(\kappa_{ij} = 1 / (1 + \tau^2 \lambda_{ij}^2)\) has a \(\mathrm{Beta}(1/2, 1/2)\) prior, the horseshoe shape that names the estimator: mass piles up at \(\kappa = 1\) hard enough to wipe out noise coefficients, and at \(\kappa = 0\) so that a large coefficient escapes shrinkage almost entirely. A Gaussian variance – the Minnesota family – shrinks everything proportionally and can do neither. There are no hyperparameters to tune, and that is real: the global scale learns the overall sparsity from the data, and the default half-Cauchy on \(\tau\) is the one the literature settled on.

Sampling uses the Makalic-Schmidt (2016) inverse-Gamma augmentation. A half-Cauchy is an inverse-Gamma mixture,

\[\lambda^2 \mid \nu \sim \mathcal{IG}(1/2, 1/\nu), \qquad \nu \sim \mathcal{IG}(1/2, 1),\]

and likewise \(\tau^2\) through an auxiliary \(\xi\), so every conditional in the hierarchy is inverse-Gamma: exact draws, no tuning, no rejection, one block per sweep. The deterministic and exogenous columns are not shrunk; they sit at a fixed loose variance, exactly where the Minnesota prior leaves them.

Variables:

deterministic_scale – Looseness of the unshrunk deterministic and exogenous rows, in units of each equation’s residual scale.

Raises:

SpecificationError – If deterministic_scale is not positive, at the first moment call, which is the model’s fit.

Note

The prior mean is zero everywhere, without a persistence hyperparameter: centring a sparsity prior away from zero would change what shrinking a coefficient means. Difference the data, or use the Minnesota family, for random-walk centring. Adaptive priors do not compose, so HorseshoePrior() + MinnesotaPrior() is refused at fit, and the Chib marginal likelihood is not offered under one; compare by predictive score instead.

Warning

Only GibbsBVAR samples the scale hierarchy. The static moments interface reports the variance at the initial scales, unit variance on the standardized lag block, which is what admissibility checks need and not the horseshoe; a model that reads only the static moments is fitting a fixed Gaussian prior under this class’s label.

See also

  • shrinkage_scales() – the posterior mean \(\tau \lambda_{ij}\) per lag coefficient, this prior’s diagnostic.

  • DirichletLaplacePrior – a global-local prior with the same interface and a Dirichlet-tied local scale.

  • NormalGammaPrior – the same interface with a tunable shape governing the tails.

  • SpikeAndSlabPrior – selection rather than shrinkage, reporting inclusion probabilities.

References

Carvalho, C. M., Polson, N. G., & Scott, J. G. (2010). The horseshoe estimator for sparse signals. Biometrika, 97(2), 465-480.

Makalic, E., & Schmidt, D. F. (2016). A simple sampler for the horseshoe estimator. IEEE Signal Processing Letters, 23(1), 179-182.

Polson, N. G., & Scott, J. G. (2010). Shrink globally, act locally: Sparse Bayesian regularization and prediction. In Bayesian Statistics 9 (pp. 501-538). Oxford University Press.

Example

Three series of which only the first has an own lag. The posterior mean local-global scale is order one where the signal is and an order of magnitude smaller on every noise coefficient:

>>> import numpy as np
>>> from cultivars.multivariate.large_dim.gibbs import GibbsBVAR
>>> rng = np.random.default_rng(0)
>>> y = np.zeros((200, 3))
>>> for t in range(1, 200):
...     y[t] = np.array([0.8, 0.0, 0.0]) * y[t - 1] + rng.standard_normal(3)
>>> res = GibbsBVAR(y, order=1, prior=HorseshoePrior()).fit(
...     n_draws=1500, n_burn=500, seed=0
... )
>>> res.prior_label
'horseshoe'
>>> scales = res.shrinkage_scales()[0]
>>> bool(scales[0, 0] > 0.8), bool(np.delete(scales, 0).max() < 0.3)
(True, True)
>>> np.abs(res.coefficients[0]).round(1)
array([[0.7, 0. , 0. ],
       [0. , 0. , 0.1],
       [0. , 0. , 0. ]])

Attributes#

deterministic_scale

Prior standard deviation of the unshrunk rows, in residual-scale units.

Methods#

coefficient_mean(context)

Zero everywhere: sparsity's claim is that coefficients are zero.

coefficient_variance(context)

The variance at the initial scales -- finite, so shrunk-admissible.

dummy_observations(context)

Artificial rows this prior contributes, for restrictions on sums.