HorseshoePrior._draw_scales#

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

Makalic-Schmidt: four inverse-Gamma blocks, all exact.

With \(m = kp \cdot k\) lag coefficients \(\beta_{ij}\) on the standardized scale, the sweep is

\[\begin{split}\lambda_{ij}^2 \mid \cdot &\sim \mathcal{IG}\bigl(1,\; 1/\nu_{ij} + \beta_{ij}^2 / (2\tau^2)\bigr), \\ \nu_{ij} \mid \cdot &\sim \mathcal{IG}\bigl(1,\; 1 + 1/\lambda_{ij}^2\bigr), \\ \tau^2 \mid \cdot &\sim \mathcal{IG}\Bigl(\tfrac{m + 1}{2},\; 1/\xi + \tfrac{1}{2} \sum_{ij} \beta_{ij}^2 / \lambda_{ij}^2\Bigr), \\ \xi \mid \cdot &\sim \mathcal{IG}\bigl(1,\; 1 + 1/\tau^2\bigr),\end{split}\]

each inverse-Gamma drawn as rate over a Gamma variate. The local scales are drawn before the global one within a sweep, and the deterministic and exogenous rows of standardized are not read.

Parameters:

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

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

  • rng (Generator) – Random generator.

  • scales (dict[str, NDArray[double]]) – The current latent state, as returned by _initial_scales() or the previous sweep.

Return type:

dict[str, NDArray[double]]

Returns:

The refreshed state under the same four names.

Example

Holding one large and three near-zero standardized coefficients fixed and sweeping the scales alone, the averaged local-global scale separates them by an order of magnitude:

>>> 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 = HorseshoePrior()
>>> rng = np.random.default_rng(0)
>>> beta = np.array([[0.0, 0.0], [3.0, 0.05], [-0.05, 0.02]])
>>> state = prior._initial_scales(ctx, prior._unit_ratio(ctx))
>>> total = np.zeros((3, 2))
>>> for _ in range(2000):
...     state = prior._draw_scales(beta, ctx, rng, state)
...     total += prior._tracked(state, ctx)
>>> (total / 2000).round(1)
array([[0. , 0. ],
       [4.2, 0.3],
       [0.3, 0.2]])