NormalGammaPrior._draw_scales#

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

GIG conditionals for the variances, Gamma for the global rate.

With \(m\) lag coefficients \(\beta_{ij}\) on the standardized scale and \(\mathrm{GIG}(p, a, b)\) the density \(x^{p-1} \exp\{-(a x + b / x) / 2\}\), the sweep is

\[\begin{split}\psi_{ij} \mid \cdot &\sim \mathrm{GIG}\bigl(\lambda - \tfrac{1}{2},\; 2\gamma,\; \beta_{ij}^2\bigr), \\ \gamma \mid \cdot &\sim \mathrm{Gamma}\Bigl(2 + m\lambda,\; 1 + \sum_{ij} \psi_{ij}\Bigr),\end{split}\]

the second in shape-rate form, the variances drawn under the incoming rate and the rate then refreshed from the new variances. Squared deviations are floored at \(10^{-11}\): at an exact zero the variance conditional is improper for \(\lambda < 1/2\), a known boundary of the hierarchy, and the floor, three parts per million on the standardized scale, restores it without moving anything statistically visible. 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; 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 refreshed state under the same two names.

Example

Holding one large and three near-zero standardized coefficients fixed and sweeping the scales alone, the averaged local 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 = NormalGammaPrior()
>>> 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. ],
       [1.8, 0.2],
       [0.2, 0.1]])