SpikeAndSlabPrior._initial_scales#

SpikeAndSlabPrior._initial_scales(context, reference)[source]#

Anchor both widths to the reference scale; start everything in.

The semiautomatic default: reference holds each coefficient’s least-squares standard error on the raw scale, and dividing by the unit ratio puts it on the standardized scale the indicators see. Spike and slab are spike and slab times that, floored at \(10^{-8}\) so a degenerate column cannot produce a zero width. Every indicator starts in the slab, so the first coefficient draw is essentially least squares, and the first sweep decides which coefficients leave.

Parameters:

  • context (_PriorContext) – The sample description; k_endog, order, scales, and lag_offset are read.

  • reference (NDArray[double]) – (width, k) reference scales, least-squares standard errors when the sampler can compute them and the unit ratio otherwise.

Return type:

dict[str, NDArray[double]]

Returns:

{"tau0", "tau1", "gamma", "prob"}, each of shape (k * p, k): the two widths, the indicators (all one), and the conditional inclusion probabilities (all inclusion).

Raises:

SpecificationError – If a hyperparameter is out of range.

Example

With a standard error of 0.05 everywhere, the spike is a two-hundredth and the slab a half:

>>> import numpy as np
>>> from cultivars.bayes.priors import PriorContext
>>> ctx = PriorContext(
...     k_endog=2, order=1, scales=np.ones(2), presample_mean=np.zeros(2)
... )
>>> state = SpikeAndSlabPrior()._initial_scales(ctx, np.full((3, 2), 0.05))
>>> state["tau0"].round(6)[0, 0], state["tau1"][0, 0], state["gamma"].sum()
(np.float64(0.005), np.float64(0.5), np.float64(4.0))