DirichletLaplacePrior._draw_scales#

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

The exact conditionals of Bhattacharya et al. (2015).

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}^{-1} \mid \cdot &\sim \mathrm{InvGaussian}\bigl(\phi_{ij} \tau / |\beta_{ij}|,\; 1\bigr), \\ \tau \mid \cdot &\sim \mathrm{GIG}\Bigl(m(a - 1),\; 1,\; 2 \sum_{ij} |\beta_{ij}| / \phi_{ij}\Bigr), \\ T_{ij} \mid \cdot &\sim \mathrm{GIG}(a - 1,\; 1,\; 2|\beta_{ij}|), \qquad \phi_{ij} = T_{ij} \Big/ \sum_{kl} T_{kl},\end{split}\]

the last line being the simplex conditional through independent GIG draws normalized to sum to one. The mixing scales are drawn under the incoming \(\phi\) and \(\tau\), then the global magnitude, then the simplex. Three guards keep the conditionals finite: \(|\beta_{ij}|\) is floored at \(10^{-10}\), the inverse-Gaussian mean is capped at \(10^8\), and the inverse-Gaussian draw is floored at \(10^{-300}\) before inversion. 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 three names.

Example

Holding one large and three near-zero standardized coefficients fixed and sweeping the scales alone, the averaged simplex puts most of its budget on the large one and the averaged local-global scale separates it 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 = DirichletLaplacePrior()
>>> 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, simplex = np.zeros((3, 2)), np.zeros((2, 2))
>>> for _ in range(2000):
...     state = prior._draw_scales(beta, ctx, rng, state)
...     total += prior._tracked(state, ctx)
...     simplex += state["phi"]
>>> (total / 2000).round(1)
array([[0. , 0. ],
       [3.8, 0.4],
       [0.4, 0.2]])
>>> bool(simplex[0, 0] / 2000 > 0.7), (simplex / 2000).sum().round(6)
(True, np.float64(1.0))