Utilities

Utilities#

Bessel Functions #

Log modified Bessel function of the second kind.

One centered whole-line Gauss–Legendre kernel (S10) implements both log_kv and log_kv_moments(). Geometry is frozen under stop_gradient; autodiff of log_kv yields cumulants of that discrete measure. No regimes, no lax.cond, no custom_jvp, no finite differences.

backend=’jax’ (default): JAX arrays, JIT-able, differentiable. backend=’cpu’ : the same sums in NumPy.

class normix.utils.bessel.BesselMoments(log_k, u0, mean, cov, d_arg, d2_arg, d2_order_arg)[source]#

Frozen quadrature moments of \((x, e^{-x}-1, e^{x}-1)\).

The discrete exponential family lives on the centered coordinate \(x = u - u_0\) with \(u_0 = \operatorname{asinh}(\nu/z)\). Geometry is constant under AD (see log_kv_moments()); mean and cov are the mean and centered Gram of \(s = (x, \operatorname{expm1}(-x), \operatorname{expm1}(x))\). cov is PSD to rounding (relative-ε entries). Argument-jet views are projections of \(w = 2\sinh(u_0 + x/2)\sinh(x/2)\), not reconstructions from cov.

Parameters:

log_k (jax.Array) – \(\log K_\nu(z)\), broadcast shape of (v, z).
u0 (jax.Array) – Mode \(\operatorname{asinh}(\nu/z)\).
mean (jax.Array) – \(E[s]\), shape (..., 3).
cov (jax.Array) – \(\mathrm{Cov}[s]\), shape (..., 3, 3).
d_arg (jax.Array) – \(\partial_z \log K_\nu = -(\cosh u_0 + E[w])\).
d2_arg (jax.Array) – \(\partial_{zz} \log K_\nu = \mathrm{Var}(w)\).
d2_order_arg (jax.Array) – \(\partial_{\nu z} \log K_\nu = -\mathrm{Cov}(x, w)\).

log_k: Array#

u0: Array#

mean: Array#

cov: Array#

d_arg: Array#

d2_arg: Array#

d2_order_arg: Array#

property d_order: Array#: \(\partial_\nu \log K_\nu = E[u] = u_0 + E[x]\).

property d2_order: Array#: \(\partial_{\nu\nu} \log K_\nu = \mathrm{Var}(u) > 0\).

normix.utils.bessel.log_kv_moments(v, z, backend='jax')[source]#

Log \(K_\nu(z)\) and the moment bundle of the quadrature measure.

One centered whole-line Gauss–Legendre kernel; backend selects the array library (JAX or NumPy), not a different approximation. Geometry is frozen under jax.lax.stop_gradient(), so the log-sum-exp identity is exact for any fixed nodes and

\[\partial_\nu \log K = E[u], \qquad \partial_{\nu\nu} \log K = \mathrm{Var}(u) > 0.\]

Stated numerical domain: \(z \ge 10^{-10}\) (GIG never asks below GIG_DEGEN_THRESHOLD), \(|\nu| \le 300\) tested. The point \((\nu,z)=(0, 10^{-300})\) is unresolved at 192 nodes.

Parameters:

v (scalar or array) – Order (any real; \(\nu=0\) is an ordinary point).
z (scalar or array) – Argument (must be \(> 0\)).
backend ('jax' or 'cpu', optional) – Array backend. 'jax' is JIT-able and differentiable; 'cpu' is the same sums in NumPy.

Returns:

log_k, u0, mean (shape (..., 3)), cov (shape (..., 3, 3)). Jet views: d_order, d_arg, d2_order, d2_order_arg, d2_arg. Argument jets (d_arg, d2_arg, d2_order_arg) are stored projections of \(w=2\sinh(u_0+x/2)\sinh(x/2)\), not affine images of cov.

Return type:

BesselMoments

Notes

jax.grad(lambda z: log_kv(v, z))(z) equals d_arg; jax.grad with respect to v equals d_order. Both are derivatives of the frozen quadrature, not finite differences. GIG uses this bundle so \(\eta\) and the Fisher come from one pass.

normix.utils.bessel.log_kv(v, z, backend='jax')[source]#

\(\log K_v(z)\) — log modified Bessel function of the second kind.

Parameters:

v (scalar or array) – Order (any real; \(K_v = K_{-v}\); \(\nu=0\) is ordinary).
z (scalar or array) – Argument (must be \(> 0\)). Stated domain \(z \ge 10^{-10}\).
backend (str, optional) – 'jax' (default) or 'cpu'. Both use the same quadrature; 'jax' is JIT-able and differentiable.

Returns:

Same broadcast shape as (v, z).

Return type:

jax.Array

Notes

A 192-point Gauss–Legendre sum of the whole-line integrand for \(K_\nu\). jax.grad with respect to z equals log_kv_moments() d_arg; with respect to v it equals d_order. Neither is a finite difference.

Examples

Evaluate at a single point (JAX backend, JIT-able). Enable float64 before importing normix:

>>> import jax
>>> jax.config.update("jax_enable_x64", True)
>>> import jax.numpy as jnp
>>> from normix import log_kv
>>> float(log_kv(v=0.5, z=1.0))
-0.774...

CPU backend (same kernel in NumPy):

>>> float(log_kv(v=0.5, z=1.0, backend='cpu'))
-0.774...

Symmetry \(K_v(z) = K_{-v}(z)\):

>>> abs(float(log_kv(0.5, 2.0)) - float(log_kv(-0.5, 2.0))) < 1e-10
True

Differentiable via JAX:

>>> import jax
>>> dlogkv_dz = jax.grad(lambda z: log_kv(0.5, z))(jnp.array(1.0))
>>> float(dlogkv_dz) < 0   # K_v decreases with z
True

Constants #

Shared numerical constants for normix.

normix.utils.constants.LOG_EPS: float = 1e-30#: Floor for log-space clamping in JAX (jnp.maximum(x, LOG_EPS))

normix.utils.constants.TINY: float = 1e-300#: Floor for numpy-side values (avoids exact zero in log)

normix.utils.constants.GIG_DEGEN_THRESHOLD: float = 1e-10#: √(ab) cap on the O(z²) remainder of the small-z GIG Gamma / InverseGamma log-partition. The limit is used iff z < this AND |p| log(2/z) exceeds BESSEL_QUAD_LOG_DROP; the form is by sign(p). Off-Θ returns +∞.

normix.utils.constants.HESSIAN_DAMPING: float = 1e-06#: Relative Tikhonov coefficient; Newton uses H + λ (tr H / n) I. Finance QP hess_reg reuses this value as an absolute ridge.

normix.utils.constants.THETA_FLOOR: float = -1e-08#: Floor for GIG θ₂, θ₃ during warm-start initialisation. Also the strict active-set gap for the Newton iteration: a coordinate within -THETA_FLOOR of a finite bound is eligible to be dropped from the KKT residual.

normix.utils.constants.KKT_NEAR_GAP: float = 0.01#: Wider gap used only to decide whether from_expectation may return a model. A GIG coordinate within this distance of θ = 0, with the gradient pointing out of the feasible set, is the bound multiplier (O(10⁻³)), not an uninverted free coordinate. The iteration itself still stops on the tighter THETA_FLOOR gap.

normix.utils.constants.BREGMAN_INVERT_ATOL: float = 1e-05#: Free residual at which η counts as inverted for from_expectation. Trust-exact status 2 stalls near 10⁻⁹ with success=False; a 20-step Newton budget on a gamma-limit GIG leaves a free residual ~10⁻⁶ and a multiplier ~10⁻³. Residuals of O(1) — the φ-gradient false stop — stay above this floor.

normix.utils.constants.GIG_THETA_PERTURB: float = 0.0001#: Near-zero perturbation for θ₂, θ₃ in GIG multi-start initialisation

normix.utils.constants.GIG_CLAMP_LO: float = 1e-06#: Clamp bounds for GIG parameters a, b in the GH M-step

normix.utils.constants.GIG_P_MAX: float = 50.0#: Maximum allowed |p| from GIG solver before falling back

normix.utils.constants.B_POST_FLOOR: float = 1e-06#: Floor for the posterior GIG scale b_post = b + (x-μ)ᵀΣ⁻¹(x-μ) in the EM E-step. Bounds the conditional inverse moment E[1/Y|x] for observations near the mode (where (x-μ)ᵀΣ⁻¹(x-μ) → 0). Only binds for VG, whose prior b = 0; for GH / NIG / NInvG the prior b > 0 already keeps b_post above this floor.

normix.utils.constants.ALPHA_MOMENT_MARGIN: float = 0.1#: Lower floor on the (α−1) denominator of the VG/NInvG prior reconstruction moment β/(α−1) in compute_eta_from_model. The Gamma inverse moment E[1/Y] (VG) and the InverseGamma forward moment E[Y] (NInvG) both equal β/(α−1) and diverge to +∞ as α↓1 (and the bare closed form turns negative for α<1). Flooring the denominator keeps the reconstructed moment finite, positive, and continuous, while leaving the well-conditioned α>1+margin regime exact. Only this single moment is regularized — E[log Y] and the other finite moment keep their exact closed forms.

normix.utils.constants.ALPHA_MIN_MARGIN: float = 0.1#: Margin ε used by the opt-in VG α-lower-bound sentinels in fit(alpha_min=…). The VG marginal density is unbounded at x=μ for α ≤ d/2, and E[1/Y|x] is unbounded for α ≤ d/2 + 1; the floored b_post keeps EM finite but cannot stop it parking near that spike. When alpha_min=’density’ the M-step clamps the Gamma shape at d/2 + ε (density bounded); alpha_min=’inverse_moment’ clamps at d/2 + 1 + ε (E[1/Y|x] also bounded). Opt-in only — default is no clamp.

normix.utils.constants.SIGMA_REG: float = 1e-08#: Regularisation added to Σ in the M-step Cholesky factorisation

normix.utils.constants.MOMENT_CANCEL_ATOL: float = 1e-06#: When ε ‖μ‖² |E[1/Y]| exceeds this, compute_eta_from_model stores the location moments about μ. A float64 uncentered second moment has already lost Σ, and no M-step arithmetic can recover it. 1e-6 matches the translation check; below it the raw moments still round-trip.

normix.utils.constants.SAFE_DENOMINATOR: float = 1e-10#: Magnitude floor for D = 1 − E[1/Y]·E[Y] in the M-step (applied as −max(|D|, floor))

normix.utils.constants.D_FLOOR: float = 1e-08#: Positivity floor for the diagonal D in the factor-analysis M-step (Σ = F Fᵀ + diag(D)). See docs/theory/factor_analysis.md.

normix.utils.constants.SIGMA_INIT_REG: float = 0.0001#: Regularisation added to empirical Σ during moment-based initialisation

normix.utils.constants.RENYI_TAYLOR_EPS: float = 1e-06#: Half-width of the Taylor window for Rényi entropy about α = 1 (H_α = H − ½ V_H (α−1) + O((α−1)²)); keeps jax.grad(renyi) defined at α = 1

normix.utils.constants.BESSEL_QUAD_NODES: int = 192#: Gauss–Legendre nodes per whole-line integral (two equal panels at the mode). Must be even. Sweep vs mpmath (n ∈ {64,96,128,192,256}): 192/40/2 is the unique setting that holds the 1e-12 scaled-error contract on the frozen table. 128/40/2 misses ∂_z at (0, 10^{-10}) (3.8e-11).

normix.utils.constants.BESSEL_QUAD_LOG_DROP: float = 40.0#: Target log-density drop at the window edge for every tilt |k| ≤ T.

normix.utils.constants.BESSEL_WINDOW_TILT: int = 2#: Largest |k| for which the window is required to capture the tilted tail.

normix.utils.constants.BESSEL_WINDOW_ITERS: int = 40#: Fixed bisection steps used to solve each window edge (jit/vmap-safe).

normix.utils.constants.BESSEL_WINDOW_HI_MAX: float = 800.0#: Absolute cap on |x| during window search (covers z down to ~1e-300).

normix.utils.constants.BESSEL_PANEL_FLOOR: float = 1e-12#: Degeneracy floor on each GL panel half-width (mode-centered).

normix.utils.constants.TORSION_SPECTRAL_FLOOR: float = 1e-12#: Relative spectral floor for Meucci torsion diagonalizations. Before correlation scaling, diag(H) is floored at TORSION_SPECTRAL_FLOOR·max(diag H); eigenvalues more negative than −TORSION_SPECTRAL_FLOOR·‖S‖_∞ mark the decomposition invalid (ENB → NaN), while roundoff-scale negatives are projected to zero. Relative so daily vs annualized returns share one value.

normix.utils.constants.MEAN_RISK_REACH_ATOL: float = 1e-08#: Infinity-norm tolerance on the constraint residual Mᵀw − c. A reduced target inside this gap is a portfolio; a larger gap is not. Attainable residuals sit at roundoff (~1e-16). The unattainable gaps in the symmetric (γ̃ = 0.1) and rank-1 (μ̃ = 0.05) examples are 0.1 and 0.03. The same absolute gap decides whether expected return is constant along a rank-2 reachable line.

Incomplete Gamma #

Pure-JAX inverse of the regularised incomplete gamma function.

JAX exposes jax.scipy.special.gammainc() (the regularised lower incomplete gamma \(P(a, x) = \gamma(a, x)/\Gamma(a)\)) but does not ship its inverse. gammaincinv() solves \(P(a, x) = q\) by Newton iteration with a Wilson–Hilferty starting guess — fully JIT- and vmap-compatible. This is the JAX analogue of scipy.special.gammaincinv().

normix.utils.gammainc.gammaincinv(a, q, max_iter=20, tol=1e-12)[source]#

Solve \(P(a, x) = q\) for \(x\) by Newton iteration.

JAX equivalent of scipy.special.gammaincinv().

Parameters:

a (Array) – Shape parameter, \(a > 0\).
q (Array) – Probability, \(q \in (0, 1)\).
max_iter (int) – Maximum Newton iterations (typically converges in <10 steps).
tol (float) – Absolute residual tolerance on \(|P(a, x) - q|\).

Return type:

Array

Notes

Starts from the Wilson–Hilferty cube-root normal approximation and iterates \(x \leftarrow x - (P(a, x) - q) / p(a, x)\), where \(p(a, x) = x^{a-1} e^{-x} / \Gamma(a)\) is the Gamma density. The density is evaluated in log space to avoid overflow for large \(a\).

PINV Sampling #

Generic RVS utilities for univariate distributions.

build_pinv_table() builds a quantile table from any univariate log-kernel in pure JAX (trapezoidal CDF on a \(w\)-grid). Distributions supply log_kernel(w) from their own log_prob (plus a Jacobian when working in \(w = \log x\)). QuantileTable holds the grids for amortised cdf / ppf / rvs; rvs_pinv() samples via inverse lookup.

class normix.utils.rvs.QuantileTable(u_grid, x_grid)[source]#

Frozen PINV quantile table (a pytree — jit/vmap/scan-safe).

Built once via build_pinv_table() (or a distribution’s quantile_table()); reuse across repeated cdf / ppf / rvs calls so the 4000-point grid is not rebuilt each time.

Parameters:

u_grid (Array)
x_grid (Array)

u_grid: Array#

x_grid: Array#

cdf(x)[source]#

CDF lookup \(F(x)\) via linear interpolation on the table.

Parameters:: x (Array)
Return type:: Array

ppf(q)[source]#

Quantile lookup \(F^{-1}(q)\) via linear interpolation.

Parameters:: q (Array)
Return type:: Array

rvs(n, seed=42)[source]#

Inverse-CDF sample of size n from this table.

Parameters:

n (int)
seed (int)

Return type:

Array

normix.utils.rvs.build_pinv_table(log_kernel, mode, *, x_of_w=None, n_grid=4000, tail_eps=1e-14)[source]#

Build a PINV quantile table in pure JAX.

Parameters:

log_kernel (Callable[[Array], Array]) – Callable w -> log f(w) for a univariate density on the internal \(w\)-axis. For support \((0, \infty)\) with \(w = \log x\), pass log_kernel(w) = log_prob(exp(w)) + w. For support \(\mathbb{R}\) with \(w = x\), pass log_kernel(w) = log_prob(w).
mode (Array) – Mode of the density on the \(w\)-axis (starting point for tail bisection).
x_of_w (Callable[[Array], Array] | None) – Map internal \(w\) to the observation axis (default identity). Use jnp.exp when \(w = \log x\).
n_grid (int) – Number of grid points for the trapezoidal CDF.
tail_eps (float) – Tail mass below which bisection stops.

Returns:

JAX arrays of shape (n_grid,) — trapezoidal CDF values and corresponding \(x\) values for jnp.interp in ppf / cdf.

Return type:

u_grid, x_grid

normix.utils.rvs.rvs_pinv(key, u_grid, x_grid, n)[source]#

Sample n observations via numerical inverse CDF.

Parameters:

key (Array) – JAX PRNG key.
u_grid (Array) – Arrays returned by build_pinv_table().
x_grid (Array) – Arrays returned by build_pinv_table().
n (int) – Sample size.

Returns:

Array of shape (n,).

Return type:

samples

Plotting #

Requires the plotting extra (uv sync --extra plotting); used by every executable tutorial for consistent figure styling.

Plotting utilities for normix notebooks.

Requires the plotting extra: pip install normix[plotting].

normix.utils.plotting.set_theme(*, scale=1.0)[source]#

Apply the normix Matplotlib theme (Kami-derived visual tokens).

Shared colour/typography tokens so tutorial figures match the documentation site aesthetic.

Parameters:: scale (float)
Return type:: None

normix.utils.plotting.style_axes(axes, *, grid_axis='y', legend=True)[source]#

Apply final styling to one or more Matplotlib axes.

Parameters:

axes (Axes | list[Axes] | tuple[Axes, ...])
grid_axis (str)
legend (bool)

Return type:

Axes | list[Axes]

normix.utils.plotting.savefig(fig, path, **kwargs)[source]#

Save a figure with normix theme defaults.

Parameters:

fig (Figure)
path (str)

Return type:

str

normix.utils.plotting.plot_pdf_cdf_comparison(configs, x, xlabel='x', title='')[source]#

Plot PDF and CDF comparing normix vs scipy for multiple distributions.

Parameters:

configs (list of dicts with keys {label, dist, scipy})
x (1-D evaluation grid)
xlabel (str)
title (str)

Return type:

Figure

normix.utils.plotting.plot_sample_histograms(configs, ncols=2)[source]#

Grid of histograms with theoretical PDF overlay.

Parameters:

configs (list of dicts with keys {label, dist, samples, x_plot (optional)})
ncols (int)

Return type:

Figure

normix.utils.plotting.plot_mle_fit(data, fit_results, xlabel='x', title='Maximum Likelihood Estimation')[source]#

Histogram of data with multiple PDF overlays for MLE comparison.

Parameters:

data (1-D sample array)
fit_results (list of dict) – Each dict has keys label, dist, ls and color.
xlabel (str)
title (str)

Return type:

Figure

normix.utils.plotting.plot_joint_1d(joint_dist, n_samples=5000, seed=42, title='')[source]#

Visualize a 1-D joint distribution f(x, y).

Left: scatter X vs Y. Right: marginal histogram of X.

Parameters:

n_samples (int)
seed (int)
title (str)

Return type:

Figure

normix.utils.plotting.plot_marginal_2d(marginal_dist, n_samples=5000, seed=42, title='')[source]#

Visualize a 2-D marginal distribution: scatter + marginal histograms.

Parameters:

n_samples (int)
seed (int)
title (str)

Return type:

Figure

normix.utils.plotting.plot_em_convergence(log_likelihoods, title='EM Convergence', true_ll=None)[source]#

Plot EM log-likelihood convergence curve.

Parameters:

log_likelihoods (List[float])
title (str)
true_ll (float | None)

Return type:

Figure

Validation #

Moment-validation and parameter-printing helpers used in notebooks and tutorials.

Moment validation and parameter printing utilities for normix notebooks.

normix.utils.validation.validate_moments(dist, n_samples=20000, seed=42, is_joint=True)[source]#

Validate E[X] and E[Y] by comparing samples vs analytical values.

Uses dist.mean() for theoretical E[X] and subordinator().mean() for E[Y].

Parameters:

n_samples (int)
seed (int)
is_joint (bool)

Return type:

Dict[str, Any]

normix.utils.validation.print_moment_validation(results, title='')[source]#

Print moment validation results.

Parameters:

results (Dict[str, Any])
title (str)

Return type:

None

normix.utils.validation.print_exp_family_params(dist, label='')[source]#

Print natural and expectation parameters.

Parameters:: label (str)
Return type:: None

Utilities

Contents

Utilities#

Bessel Functions#

Constants#

Incomplete Gamma#

PINV Sampling#

Plotting#

Validation#

Bessel Functions #

Constants #

Incomplete Gamma #

PINV Sampling #

Plotting #

Validation #