SpikeAndSlabPrior._draw_scales#

SpikeAndSlabPrior._draw_scales(standardized, context, rng, scales)[source]#

Exact Bernoulli conditionals for the indicators.

For each lag coefficient the log odds of inclusion are

\[\log\frac{p}{1 - p} + \log\frac{\tau_0}{\tau_1} + \frac{\beta^2}{2}\Bigl(\frac{1}{\tau_0^2} - \frac{1}{\tau_1^2}\Bigr),\]

clipped to \(\pm 700\) before the logistic so a coefficient far outside the spike cannot overflow, and the indicator is a Bernoulli draw at that probability. The widths do not change across sweeps; only gamma and prob are refreshed. The probability is kept alongside the draw because it, not the binary indicator, is what _tracked() averages.

Parameters:

  • standardized (NDArray[double]) – (width, k) current coefficients divided by the unit ratio.

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

  • rng (Generator) – Random generator.

  • scales (dict[str, NDArray[double]]) – The current latent state.

Return type:

dict[str, NDArray[double]]

Returns:

The state with gamma and prob refreshed.

Example

A coefficient twelve spike widths from zero is in with certainty; one well inside the spike is in with roughly the prior odds discounted by the width ratio:

>>> 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)
... )
>>> prior = SpikeAndSlabPrior()
>>> state = prior._initial_scales(ctx, np.full((3, 2), 0.05))
>>> beta = np.array([[0.0, 0.0], [0.6, 0.001], [-0.03, 0.001]])
>>> state = prior._draw_scales(beta, ctx, np.random.default_rng(0), state)
>>> state["prob"].round(3)
array([[1.  , 0.01],
       [1.  , 0.01]])
>>> state["gamma"]
array([[1., 0.],
       [1., 0.]])