sys_mapping.model_selection

Likelihood ratio test (LRT) for comparing nested contamination models.

likelihood_ratio_test() computes \(\lambda_{\rm LR} = 2(\ln\mathcal{L}_{\rm alt} - \ln\mathcal{L}_{\rm null})\) and returns the p-value against a \(\chi^2(r)\) distribution where \(r\) is the difference in free parameters.

Typical use: test additive (null) vs. combined (alternative) to decide whether multiplicative contamination is significant.

Key paper: Berlfein et al. 2024, Eq. 19 — see also Methods Reference.

Likelihood ratio test for model selection (Eq. 19, Berlfein et al. 2024).

The likelihood ratio statistic:: λ_LR = 2 [ln L(Θ̂) - ln L(Θ̂_0)]

is asymptotically χ²(r) distributed under the null hypothesis, where r = dim(Θ̂) - dim(Θ̂_0) is the number of additional free parameters.

class sys_mapping.model_selection.LikelihoodRatioResult(lambda_lr: 'float', n_dof: 'int', p_value: 'float', reject_null: 'bool', null_model: 'str', alt_model: 'str')[source]

Bases: object

Parameters:

lambda_lr (float)
n_dof (int)
p_value (float)
reject_null (bool)
null_model (str)
alt_model (str)

lambda_lr: float

n_dof: int

p_value: float

reject_null: bool

null_model: str

alt_model: str

sys_mapping.model_selection.likelihood_ratio_test(delta_g_obs, delta_t, theta_null, theta_alt, null_model, alt_model, use_skewed=False, significance=0.05)[source]

Perform a likelihood ratio test comparing null vs. alternative model (Eq. 19).

Computes λ_LR = 2 [ln L(θ̂_alt) − ln L(θ̂_null)], which under the null hypothesis is asymptotically χ²(r) where r = dim(alt) − dim(null) is the extra degrees of freedom.

Note

Each call creates new JIT-compiled likelihood functions internally. For repeated calls with the same data, precompute the log-likelihoods using make_log_likelihood() to avoid recompilation.

Parameters:

delta_g_obs ((n_pix,) observed overdensity)
delta_t ((n_sys, n_pix) template values)
theta_null (flat parameter vector for the null model (from get_mle_params()))
theta_alt (flat parameter vector for the alternative model)
null_model (str 'additive' or 'multiplicative' (must be nested in alt_model))
alt_model (str 'combined' (must have more free parameters than null_model))
use_skewed (bool)
significance (float p-value threshold for rejecting the null; default 0.05)

Returns:

LikelihoodRatioResult – lambda_lr : test statistic (≥ 0 only when both theta are at their MLEs) n_dof : degrees of freedom r = n_free_params(alt) − n_free_params(null) p_value : Pr(χ²(r) > λ_LR) under the null reject_null: True when p_value < significance null_model, alt_model: model names
Performance
———–
Measured on CPU (n_pix=10_000, n_sys=3) (*~481 ms/call (includes JIT*)
compilation of two likelihood functions, ~227 ms each on first call).
In a pipeline where likelihoods are already compiled, evaluation is ~1 ms.
Precision
———
p-value is in [0, 1] by construction.
lambda_lr = 0 when theta_alt uses b=0 (equivalent to the null).
For large n_pix, lambda_lr is χ²-distributed under the null asymptotically.

Return type:

LikelihoodRatioResult

Examples

>>> import numpy as np
>>> from sys_mapping import likelihood_ratio_test, pack_params
>>> rng = np.random.default_rng(0)
>>> n_pix, n_sys = 5000, 3
>>> delta_g = rng.standard_normal(n_pix) * 0.1
>>> delta_t = rng.standard_normal((n_sys, n_pix))
>>> a_hat = np.zeros(n_sys)
>>> b_hat = np.zeros(n_sys)
>>> sigma  = float(delta_g.std())
>>> theta_add  = pack_params(a_hat, None,  sigma, model="additive")
>>> theta_comb = pack_params(a_hat, b_hat, sigma, model="combined")
>>> result = likelihood_ratio_test(delta_g, delta_t,
...                                theta_add, theta_comb,
...                                "additive", "combined")
>>> result.n_dof   # combined has n_sys extra b parameters
3
>>> 0.0 <= result.p_value <= 1.0
True

class sys_mapping.model_selection.ForwardSelectionRound(round_num, added_index, p_value)[source]

Bases: object

Diagnostics for one round of greedy forward template selection.

Parameters:

round_num (int)
added_index (int)
p_value (float)

round_num: int

added_index: int

p_value: float

class sys_mapping.model_selection.GreedyForwardSelectionResult(selected_indices, rounds, p_threshold, n_initial)[source]

Bases: object

Result of greedy_forward_select().

Parameters:

selected_indices (list[int])
rounds (list[ForwardSelectionRound])
p_threshold (float)
n_initial (int)

selected_indices

Indices into the input delta_t of accepted templates, in selection order.

Type:: list[int]

rounds

Per-round diagnostics: which template was added and at what p-value.

Type:: list[ForwardSelectionRound]

p_threshold

The LRT gate used.

Type:: float

n_initial

Number of candidate templates passed in.

Type:: int

selected_indices: list[int]

rounds: list[ForwardSelectionRound]

p_threshold: float

n_initial: int

sys_mapping.model_selection.greedy_forward_select(delta_g, delta_t, p_threshold=0.05, use_skewed=False)[source]

Greedy LRT-based forward template selection.

Starting from an empty model, iteratively adds the template whose inclusion yields the smallest χ²(1) LRT p-value, stopping when no remaining candidate achieves p_value < p_threshold. OLS is used to compute the MLE at each step, avoiding the cost of full MCMC inference during selection.

Parameters:

delta_g ((n_pix,) observed galaxy overdensity.)
delta_t ((n_cand, n_pix) candidate template pixel values.)
p_threshold (float) – Significance gate; a candidate is added only if its p-value is strictly below this threshold. Default 0.05.
use_skewed (bool) – Whether to include a skewness parameter in the likelihood. Default False.

Returns:

.selected_indices lists the row indices into delta_t for the accepted subset, in the order they were added.

Return type:

GreedyForwardSelectionResult

Notes

JIT-compiled likelihood functions are cached at module level the first time each (n_sys, use_skewed) combination is encountered, so subsequent calls with the same sizes are fast.

Examples

>>> import numpy as np
>>> from sys_mapping import greedy_forward_select
>>> rng = np.random.default_rng(0)
>>> n_pix = 2000
>>> delta_t = rng.standard_normal((5, n_pix))
>>> # Pure noise — no template should be significant at very tight threshold
>>> result = greedy_forward_select(rng.standard_normal(n_pix), delta_t,
...                                p_threshold=1e-6)
>>> len(result.selected_indices) == 0
True

class sys_mapping.model_selection.SnrPreselectionResult(selected_indices, snr_values, method, snr_min, n_top, n_initial)[source]

Bases: object

Result of snr_preselect().

Parameters:

selected_indices (list[int])
snr_values (ndarray)
method (str)
snr_min (float | None)
n_top (int | None)
n_initial (int)

selected_indices

Indices into delta_t, sorted from highest to lowest SNR.

Type:: list[int]

snr_values

SNR/Δχ² for all n_initial templates (shape (n_initial,)).

Type:: np.ndarray

method

SNR estimator used ("data", "template", "peak", or "isd").

Type:: str

snr_min

Minimum SNR threshold applied, or None if not used.

Type:: float or None

n_top

Top-K cap applied, or None if not used.

Type:: int or None

n_initial

Total number of candidate templates passed in.

Type:: int

selected_indices: list[int]

snr_values: ndarray

method: str

snr_min: float | None

n_top: int | None

n_initial: int

sys_mapping.model_selection.snr_preselect(delta_g_obs, delta_t, *, method='data', snr_min=None, n_top=None, n_bins=10, poly_order=1, fracdet=None)[source]

Rank and pre-select templates by SNR of their cross-correlation with data.

Step 1 of the two-stage decontamination pipeline: compute a cheap SNR metric for every candidate template and keep only those above a threshold and/or in the top-K by SNR. Pass the selected_indices output to greedy_forward_select() (or direct MCMC) for Step 2.

For mock-based significance (the paper method), call isd_template_significance() directly and threshold its "p_values" output.

Parameters:

delta_g_obs ((n_pix,) observed galaxy overdensity.)
delta_t ((n_cand, n_pix) candidate template pixel values.)
method (str) – SNR estimator; one of "data", "template", "peak", "isd". Default "data" (Pearson cross-correlation).
snr_min (float | None) – Keep only templates with snr >= snr_min. If None, no threshold is applied.
n_top (int | None) – Keep at most the n_top highest-SNR templates. Applied after snr_min. If None, no cap is applied.
n_bins (int) – Number of equal-width bins for method="isd".
poly_order (int) – Polynomial degree for the 1D fit in method="isd".
fracdet (ndarray | None) – Per-pixel fractional coverage weights for method="isd".

Returns:

.selected_indices lists the accepted template indices in descending SNR order. .snr_values contains the SNR for all n_initial templates (useful for inspection).

Return type:

SnrPreselectionResult

Examples

>>> import numpy as np
>>> from sys_mapping import snr_preselect, greedy_forward_select
>>> rng = np.random.default_rng(0)
>>> n_pix, n_cand = 4000, 8
>>> delta_t = rng.standard_normal((n_cand, n_pix))
>>> delta_g = 0.6 * delta_t[3] + rng.standard_normal(n_pix) * 0.1
>>> pre = snr_preselect(delta_g, delta_t, method="data", n_top=4)
>>> len(pre.selected_indices)
4
>>> pre.selected_indices[0]  # template 3 should rank first
3
>>> pre.snr_values.shape
(8,)