pmcprg.diagnostics

Goodness-of-fit and inference diagnostics, independent of any specific model class: each utility takes raw NumPy arrays (and optionally a CDF callable) and returns a typed result object.

Multivariate Kolmogorov-Smirnov test

Naaman (2021)'s finite-sample extension.

mks_1samp

mks_1samp(sample: ndarray, cdf: Callable[[ndarray], float], *, alpha: float = 0.05, asymptotic: bool = False) -> MKSResult

One-sample multivariate KS test.

Test H0: sample is drawn from the distribution whose CDF is cdf.

Parameters:

Name	Type	Description	Default
sample	ndarray		required
cdf	Callable[[ndarray], float]	return a scalar in [0, 1] (the joint CDF evaluated at the point).	required
alpha	float		0.05
asymptotic	bool — use the asymptotic critical value (tighter, but	only valid for large ``N``). Default False (finite-sample Naaman bound).	False

Returns:

Type	Description
MKSResult — see :class:`MKSResult`.

mks_2samp

mks_2samp(sample_a: ndarray, sample_b: ndarray, *, alpha: float = 0.05, asymptotic: bool = False) -> MKSResult

Two-sample multivariate KS test.

Test H0: sample_a and sample_b share the same distribution.

Parameters:

Name	Type	Description	Default
sample_a	np.ndarray, shape ``(N_a, d)`` and ``(N_b, d)``.	Must have the same trailing dimension. 1-D arrays are reshaped to (N, 1) automatically.	required
sample_b	np.ndarray, shape ``(N_a, d)`` and ``(N_b, d)``.	Must have the same trailing dimension. 1-D arrays are reshaped to (N, 1) automatically.	required
alpha	float		0.05
asymptotic	float		0.05

Returns:

Type	Description
MKSResult — see :class:`MKSResult`.

mks_test

mks_test(sample: ndarray, other: ndarray | None = None, cdf: Callable[[ndarray], float] | None = None, *, alpha: float = 0.05, asymptotic: bool = False) -> MKSResult

Dispatch to :func:mks_1samp or :func:mks_2samp.

Exactly one of other / cdf must be provided.

Parameters:

Name	Type	Description	Default
sample	see :func:`mks_1samp` and :func:`mks_2samp`.		required
other	see :func:`mks_1samp` and :func:`mks_2samp`.		required
cdf	see :func:`mks_1samp` and :func:`mks_2samp`.		required
alpha	float		0.05
asymptotic	float		0.05

Returns:

Type	Description
MKSResult.

MKSResult dataclass

MKSResult(statistic: float, critical_value: float, alpha: float, reject: bool, n_samples_x: int, n_samples_y: int | None, dim: int, asymptotic: bool)

Typed result of a multivariate Kolmogorov-Smirnov test.

Fields

statistic : float — the KS statistic (max absolute deviation). critical_value : float — the critical value at significance level :attr:alpha. Computed under the finite-sample bound (default) or the asymptotic approximation (asymptotic=True). alpha : float — the significance level used. reject : bool — True if statistic > critical_value (H0 rejected at level α). n_samples_x : int — size of the (first) sample. n_samples_y : int | None — size of the second sample (None for the one-sample test). dim : int — observation dimensionality. asymptotic : bool — whether the critical value uses the asymptotic approximation.

Parametric bootstrap

Model-level calibration of any post-fit statistic — resampling whole series from the fitted model, needed because a tabulated critical value is wrong on a fitted model (estimated reference, dependent observations, per-state samples that are assigned rather than given).

parametric_bootstrap

parametric_bootstrap(model, Y, statistic: Callable[[Any, ndarray, Generator], Mapping[Hashable, float]], *, B: int = 200, seed: int = 0, max_len: int | None = None, progress_cb: Callable[[int, int, float, str], None] | None = None, simulate_fn: Callable[..., tuple] | None = None, n_jobs: int | None = None) -> dict[Hashable, BootstrapResult]

Calibrate statistic against series simulated from model.

Parameters:

Name	Type	Description	Default
model	PMCModel	The fitted model. It is both what is being tested and what generates the null — see the module docstring on what that means for the hypothesis actually being tested.	required
Y	array - like	The observations, shape (N,) or (N, d).	required
statistic	callable	statistic(model, Y, rng) -> {key: float}. Called once on the real data and once per replicate. It may consume randomness — a posterior draw, for instance — and gets a generator for that; observed and replicates draw from the same stream, so a statistic that uses one draw is compared against nulls that use one draw too.	required
B	int	Replicates. Cost is linear in N·B.	200
max_len	int	Cap on the series length. Above it a contiguous prefix is used, and replicates are simulated at that same length so calibration is untouched — what is lost is sample size, not validity. Contiguous because thinning would destroy the serial dependence the bootstrap exists to absorb.	None
progress_cb	callable	progress_cb(b, B, 0.0, "bootstrap") before each replicate.	None
simulate_fn	callable	Defaults to :func:pmcprg.pmc.simulate.simulate. Injectable for tests and for model classes that simulate differently.	None
n_jobs	int	None (default) keeps the historical loop: one generator, default_rng(seed), feeds the observed pass and then every replicate in turn — its simulation seed, then whatever statistic draws from it. What replicate b receives depends on how many numbers the replicates before it consumed (a posterior draw consumes a data-dependent amount), so that stream cannot be split across processes and stays the default, unchanged to the bit. An int switches to per-replicate streams (audit FR-14): replicate b gets its own generator, default_rng(SeedSequence(seed).spawn(B)[b]), draws its simulation seed from it and passes it to statistic. Replicates no longer depend on one another nor on the order they run in, so n_jobs=1 (in this process), k (k worker processes) and -1 (one per CPU) return identical results — but not those of n_jobs=None: same law, other numbers. The observed pass still draws from default_rng(seed) and is unchanged. With k > 1, model, statistic and simulate_fn are sent to the workers and must be picklable: top-level functions of an importable module, not lambdas or closures (TypeError otherwise). progress_cb is then called with the number of replicates finished so far: 0 first, then after each replicate (in-process) or each chunk (with workers). Cost: starting the workers takes about a second, then one simulation and one statistic per replicate. With the GUI's copula statistic (forward-backward, one ξ-weighted CvM per pair) on a two-state PMC, B = 200: N = 800, 13.3 s → 4.4 s with 4 workers, 3.4 s with 8; N = 2000, 40.7 s → 11.3 s and 7.7 s (Apple M2 Pro, 6 performance + 4 efficiency cores; full table in the CHANGELOG entry of FR-14). Warnings and log records raised in the workers are re-emitted here (:mod:pmcprg._parallel).	None

Returns:

Type	Description
dict mapping each key returned by ``statistic`` to a
class:`BootstrapResult`. Keys seen only in replicates and never in the
observed pass are dropped: there is nothing to compare them to.

BootstrapResult dataclass

BootstrapResult(observed: float, p_value: float, n_valid: int, B: int, n_series: int)

One statistic's observed value and its bootstrap p-value.

Fields

observed : float — the statistic on the real data. p_value : float — (1 + #{T_b ≥ T_obs}) / (n_valid + 1), the exactly valid Monte-Carlo p-value (Davison & Hinkley 1997, ch. 4; Phipson & Smyth 2010): it can never be 0, and the smallest reportable value is 1/(B + 1). NaN when no replicate landed. n_valid : int — replicates that produced a finite value. Report it: a p-value resting on 12 of 200 draws is not the same claim as one resting on 200, and the difference is invisible in the number itself. B : int — replicates requested. n_series : int — length of the series actually used, which may be shorter than the input (see max_len).

Smooth tests (Neyman)

Legendre-polynomial components of a fitted margin's goodness of fit — one component per kind of departure (location, dispersion, asymmetry, tails).

neyman_components

neyman_components(sample, cdf=None, *, order: int = 4) -> NeymanResult

Smooth-test components of sample against cdf.

Parameters:

Name	Type	Description	Default
sample	array-like, shape ``(n,)``	The observations. One-dimensional: the construction rests on the probability integral transform, which has no direct multivariate analogue here.	required
cdf	callable	Vectorised CDF of the reference distribution. Pass None when sample already holds the transformed values.	None
order	int	Number of components, at most 4 — beyond V₄ the polynomials answer questions no margin diagnostic asks, and each one costs multiplicity.	4

Returns:

Type	Description
NeymanResult

neyman_test

neyman_test(p_values, *, alpha: float = 0.05) -> tuple[float, int]

Combine per-component p-values into one verdict, paying multiplicity.

p_values are the bootstrap p-values of V_j², in component order. Returns (corrected_p, index) where index is the 1-based position of the component that drove it.

Bonferroni, not the omnibus Σ V_j². Both hold their level; measured over 200 datasets per alternative, Bonferroni rejects a Student-t margin at 0.338 against the omnibus's 0.180 and a skew-normal one at 0.480 against 0.290. Concentrating on one component and paying a factor of four beats spreading over four degrees of freedom, because a real margin is usually wrong in one way.

Reporting the uncorrected minimum would take the level from 0.05 to 0.175 — measured on the same datasets. The factor is not optional.

NeymanResult dataclass

NeymanResult(components: ndarray, squares: ndarray, n: int, component: str = '', index: int = 0)

Components of a smooth test, and the multiplicity-corrected verdict.

Fields

components : np.ndarray — the signed V_j, asymptotically N(0, 1) under a fully specified null; with an estimated margin the null law changes (Kallenberg & Ledwina 1997), which is why the package calibrates them by parametric_bootstrap. squares : np.ndarray — V_j², the per-component statistics. n : int — sample size the components were computed on. component : str — name of the largest component: where the margin departs, which a single p-value cannot say. index : int — its 1-based position.

Read the signed components, not only the largest. A single name is a summary and summaries mislead: heavy tails show up as V₂ negative and V₄ positive, because extra mass in the extremes comes with a hollowed-out shoulder. Taking the argmax alone would label that "dispersion" and lose half the diagnosis.

COMPONENT_NAMES module-attribute

COMPONENT_NAMES = ('location', 'dispersion', 'asymmetry', 'tail weight')

Copula-family comparison (Vuong / Clarke)

Is a copula family choice significant? Weighted Vuong (1989) and Clarke (2007) tests with a HAC variance, over all pairs of a candidate set, and over every pair of states after an ICE run.

vuong_test

vuong_test(logc_a, logc_b, weights=None, *, k_a: float = 1, k_b: float = 1, correction: str = 'none', alpha: float = 0.05, bandwidth: int | str = 'auto', n_eff: str = 'sum', names: tuple[str, str] = ('A', 'B'), ranks: bool = False, uv=None, dm_du=None, dm_dv=None, pretest_omega2: bool = False) -> ComparisonResult

Weighted Vuong (1989) test of family A against family B.

Parameters:

Name	Type	Description	Default
logc_a	pointwise log-densities ``log c_A(u_n, v_n; θ̂_A)`` and	``log c_B(u_n, v_n; θ̂_B)`` at each family's fitted parameters, shape ``(N,)``.	required
logc_b	pointwise log-densities ``log c_A(u_n, v_n; θ̂_A)`` and	``log c_B(u_n, v_n; θ̂_B)`` at each family's fitted parameters, shape ``(N,)``.	required
weights			None
k_a	float		1
k_b	float		1
correction	str		'none'
alpha	float		0.05
bandwidth	int | str	(Newey–West 1994, at least 1).	'auto'
n_eff	str		'sum'
names	tuple[str, str]		('A', 'B')
ranks	bool	the variance (default ``False``, the historical behaviour — margins treated as known, as documented in the "Estimated pseudo-observations" caveat of the module docstring). Needs ``uv``, ``dm_du``, ``dm_dv`` (see :func:`margin_correction_derivatives`); the point estimate ``m̄_w`` is unaffected, only ``σ̂²`` is. Adapted from the ``W₁, W₂`` sandwich correction of ``pmcprg.copulas._stderr.standard_errors(ranks=True)``, applied to the log-density *difference* instead of one family's score — the same kind of extra variance term from the margins' own estimation, not a term re-derived from Chen & Fan's paper directly (offline; see the module docstring's honesty note and ``test_model_selection_mc.py`` for the measured level improvement on rank pseudo-observations).	False
uv	required when ``ranks=True`` — the pseudo-observations	and the pointwise ``∂m/∂u``, ``∂m/∂v`` from :func:`margin_correction_derivatives`.	None
dm_du	required when ``ranks=True`` — the pseudo-observations	and the pointwise ``∂m/∂u``, ``∂m/∂v`` from :func:`margin_correction_derivatives`.	None
dm_dv	required when ``ranks=True`` — the pseudo-observations	and the pointwise ``∂m/∂u``, ``∂m/∂v`` from :func:`margin_correction_derivatives`.	None
pretest_omega2	run Vuong's (1989) ``H0: ω² = 0`` pre-test (see	:func:`omega2_test`; not the exact weighted-χ² construction, a documented substitute) and attach it as ``result.omega2``; when it does **not** reject, the decision is short-circuited to ``"tie"`` (method ``"omega2_pretest"``) instead of trusting the normal Z-test, following Vuong's own recommendation. Default ``False``: unchanged behaviour.	False

Returns:

Type	Description
ComparisonResult — ``Z = √n_eff (m̄_w − K/n_eff) / σ̂``, two-sided normal
p-value, decision ``"A"`` / ``"B"`` / ``"tie"``. See the module docstring
for the variance, the corrections and the caveats (margins estimated,
Chen & Fan 2006; overlapping families).

clarke_test

clarke_test(logc_a, logc_b, weights=None, *, k_a: float = 1, k_b: float = 1, correction: str = 'none', alpha: float = 0.05, bandwidth: int | str = 'auto', n_eff: str = 'sum', method: str = 'auto', names: tuple[str, str] = ('A', 'B')) -> ComparisonResult

Weighted distribution-free sign test of Clarke (2007).

H₀: median(m_n) = 0, with the corrected differences m_n − K / n_eff (Clarke 2007). Points whose corrected difference is exactly zero carry no sign and are dropped, as in the classical sign test.

Takes the same parameters as :func:vuong_test, plus method.

Parameters:

Name	Type	Description	Default
method	``"exact"`` — binomial ``B ~ Bin(N, 1/2)`` on the weighted count	``B = Σ w 1{m > 0}`` of ``N = Σ w`` trials; valid only with integer (frequency) weights, ``bandwidth=0`` and ``n_eff="sum"``. ``"normal"`` — ``Z = √n_eff (p̂ − 1/2) / σ̂``, ``p̂ = B/N``, with the HAC variance of ``1{m > 0} − 1/2`` (centred at the null value, so that ``bandwidth=0`` gives ``σ̂² = 1/4``, the binomial variance; the automatic bandwidth is chosen on the series centred at ``p̂``). ``"auto"`` — ``"exact"`` when its conditions hold, ``"normal"`` otherwise.	'auto'

Notes

The binomial law assumes independent signs: with serially dependent pairs (ICE) use the default HAC.

comparison_matrix

comparison_matrix(log_densities: Mapping[str, ndarray], weights=None, *, n_params: Mapping[str, float] | None = None, test: str = 'vuong', correction: str = 'none', alpha: float = 0.05, bandwidth: int | str = 'auto', n_eff: str = 'sum', **test_kwargs) -> ComparisonMatrix

Pairwise tests between all families of a candidate set.

Parameters:

Name	Type	Description	Default
log_densities	``{name: log c_name(u_n, v_n; θ̂_name)}``, every array on	the same ``N`` points, in the order of the output.	required
weights			None
n_params	Mapping[str, float] | None		None
test	str		'vuong'
**test_kwargs	forwarded to the test.		{}

Notes

Each unordered pair is tested once; the lower triangle is the swapped result, so the matrix is exactly antisymmetric.

confidence_set

confidence_set(log_densities: Mapping[str, ndarray], weights=None, *, n_params: Mapping[str, float] | None = None, best: str | None = None, scores: Mapping[str, float] | None = None, test: str = 'vuong', correction: str = 'none', alpha: float = 0.05, adjust: str | None = None, bandwidth: int | str = 'auto', n_eff: str = 'sum', **test_kwargs) -> ConfidenceSet

Families not significantly worse than the best one.

The reference best is, by default, the family with the highest scores — or, without scores, the highest corrected weighted log-likelihood Σ w log c − K (the ranking of AIC for correction="akaike", of BIC for "schwarz"). Every other family F is tested against it; F is excluded only when the test significantly prefers best (p ≤ alpha after adjust and a statistic of the sign favouring best). The rest are retained: they are statistically tied with the best, or better than it under the test (possible when best comes from a non-likelihood score).

adjust (None, "holm", "bonferroni") corrects the m − 1 p-values for multiplicity. This is a set of pairwise comparisons with one reference, not the model confidence set of Hansen, Lunde & Nason (2011).

ice_pair_comparisons

ice_pair_comparisons(model, Y, xi=None, *, candidates: Sequence[str] | None = None, criterion: str | None = None, test: str = 'vuong', correction: str | None = None, alpha: float = 0.05, adjust: str | None = None, bandwidth: int | str = 'auto', n_eff: str = 'sum', pairs: Sequence[tuple[int, int]] | None = None, **test_kwargs) -> dict[tuple[int, int], PairComparison]

Is the copula family ICE selected for each pair of states significant?

For every pair (i, j) this rebuilds the weighted sample the ICE M-step fits the copula c_ij on — (F_ij(y_n), F_ji(y_{n+1})) with weight ξ_n(i, j) (DerrodePieczynski_CSDA2013 Eq. 12), the CDFs clipped as in pmcprg.pmc.ice._m_step — refits every candidate with the M-step's own weighted MLE (ice._fit_copula_params), scores it with the ICE criterion (ice._SCORE_FN), and compares the candidates with :func:vuong_test or :func:clarke_test.

Parameters:

Name	Type	Description	Default
model			required
Y			required
xi		(forward–backward) of ``model`` on ``Y``, as ICE does.	None
candidates	copula SHORT_NAMEs; ``None`` uses the ICE configuration of	the model (``[ice]`` section, else the package default).	None
criterion	str | None	uses the model's ICE configuration.	None
test	str		'vuong'
correction	``None`` takes the one consistent with ``criterion``	(``mle`` → ``"none"``, ``aic`` → ``"akaike"``, ``bic`` → ``"schwarz"``, other criteria → ``"none"``).	None
alpha	see :func:`confidence_set`.		0.05
adjust	see :func:`confidence_set`.		0.05
bandwidth	see :func:`confidence_set`.		0.05
n_eff	see :func:`confidence_set`.		0.05
pairs	Sequence[tuple[int, int]] | None	weight is not negligible (ICE skips the others too).	None

Returns:

Type	Description
``{(i, j): PairComparison}``. ``tied_with_runner_up`` answers the
question "is the selected family significantly better than the next
one?". Nothing in ``model`` is modified.

Notes

selected is the criterion's choice on this sample; it equals the model's family when ICE stopped at convergence (the returned model then produced ξ), and may differ after max_iter iterations. The test ignores the estimation of the margins and of ξ (Chen & Fan 2006; see the module docstring).

ComparisonMatrix dataclass

ComparisonMatrix(names: list[str], statistic: ndarray, p_value: ndarray, decision: list[list[str]], results: dict[tuple[str, str], ComparisonResult], test: str, alpha: float, correction: str)

All pairwise comparisons of a candidate set.

statistic[a, b] is the statistic of names[a] against names[b] (antisymmetric, 0 on the diagonal); p_value is symmetric (1 on the diagonal); decision[a][b] is the preferred family's name, "tie" or "undetermined" ("—" on the diagonal). results[(A, B)] holds the :class:ComparisonResult for every ordered pair A ≠ B.

ties

ties() -> list[tuple[str, str]]

Unordered pairs that the test does not separate.

ComparisonResult dataclass

ComparisonResult(test: str, names: tuple[str, str], statistic: float, p_value: float, alpha: float, decision: str, method: str, estimate: float, penalty: float, sd: float, n_eff: float, sum_weights: float, bandwidth: int, correction: str, n_eff_kind: str, count: float = math.nan, n_trials: float = math.nan, ranks: bool = False, omega2: 'Omega2Test | None' = None)

Outcome of one pairwise test of family A against family B.

Attributes:

Name	Type	Description
test	``"vuong"`` or ``"clarke"``.
names	``(name_A, name_B)``.
statistic	standardised statistic, ``> 0`` favours ``A``. For Clarke's	exact test, (B − N/2)/√(N/4) with B the count.
p_value	two-sided p-value.
alpha	level of the decision.
decision	``"A"``, ``"B"``, ``"tie"`` or ``"undetermined"``.
method	``"normal"``, ``"exact"`` (binomial), ``"infinite"`` (one	likelihood is −∞), "degenerate" (zero variance or no informative point), "undetermined".
estimate	Vuong — weighted mean ``m̄_w`` of the *uncorrected*	differences (nat per unit weight); Clarke — weighted share of positive corrected differences.
penalty	total correction ``K`` subtracted from ``Σ w m``.
sd	long-run standard deviation ``σ̂`` (NaN when undefined).
n_eff	effective sample size (``Σw`` or Kish).
sum_weights	``Σ w`` over the points that entered the test.
bandwidth	HAC lag truncation ``L`` actually used.
correction, n_eff_kind	the options used.
count, n_trials	Clarke only — weighted count of positive differences and	weighted number of non-zero differences.
ranks	Vuong only — whether the Chen & Fan (2006) rank-margin	correction was applied (see :func:vuong_test's ranks).
omega2	Vuong only, when ``pretest_omega2=True`` — the	:class:Omega2Test of H0: ω² = 0 that gated the decision.

tie property

tie: bool

True when the test does not separate the two families.

preferred property

preferred: str | None

Name of the significantly preferred family, None otherwise.

swapped

swapped() -> 'ComparisonResult'

The same comparison read as B against A.

ConfidenceSet dataclass

ConfidenceSet(best: str, members: list[str], excluded: list[str], scores: dict[str, float], comparisons: dict[str, ComparisonResult], p_adjusted: dict[str, float], alpha: float, adjust: str | None, test: str, correction: str)

Families not significantly worse than a reference family.

Attributes:

Name	Type	Description
best	the reference family (kept by construction).
members	retained families, reference first, then by decreasing score.
excluded	families significantly worse than ``best``.
scores	the score used for the ordering (corrected weighted	log-likelihood, or the caller's scores).
comparisons	``{name: test(best, name)}`` for every other family.
p_adjusted	the p-values after ``adjust`` (equal to the raw ones when	adjust is None).

is_tie property

is_tie: bool

True when at least one other family is retained with best.

PairComparison dataclass

PairComparison(pair: tuple[int, int], model_family: str, selected: str, runner_up: str | None, criterion: str, scores: dict[str, float], log_likelihood: dict[str, float], params: dict[str, dict], n_params: dict[str, int], sum_weights: float, matrix: ComparisonMatrix | None, vs_runner_up: ComparisonResult | None, confidence_set: ConfidenceSet | None, failed: list[str] = list(), log_densities: dict[str, ndarray] = dict(), weights: ndarray | None = None)

Family comparison for one pair of states (i, j) of a fitted chain.

Attributes:

Name	Type	Description
pair	``(i, j)``.
model_family	the family the fitted model carries for the pair.
selected	the family the ICE criterion selects on this weighted	sample (what the next M-step would keep; equal to model_family at a converged run).
runner_up	second family by the criterion (``None`` with a single	usable candidate).
criterion	the ICE selection criterion used for the ranking.
scores	criterion score of every candidate (higher is better).
log_likelihood	``Σ ξ log c`` of every candidate at its fitted parameters.
params	fitted parameters of every candidate.
n_params	parameter count of every candidate.
sum_weights	``Σ_n ξ_n(i, j)``.
matrix	:class:`ComparisonMatrix` over the usable candidates.
vs_runner_up	``test(selected, runner_up)``.
confidence_set	families not significantly worse than ``selected``.
failed	candidates whose fitted copula could not be evaluated.
log_densities	``{family: log c(u_n, v_n)}`` at the fitted parameters, on	the N − 1 transitions — to run another test without refitting, e.g. clarke_test(ld[a], ld[b], weights).
weights	``ξ_n(i, j)``, shape ``(N − 1,)``.

tied_with_runner_up property

tied_with_runner_up: bool

True when the selected family is not significantly better than the runner-up (False when there is no runner-up).

Radial symmetry and exchangeability screens

Cramér-von Mises statistics on the empirical copula against its radial reflection / its transpose, calibrated by a parametric (Gaussian-surrogate) bootstrap. Unweighted only.

radial_symmetry_test

radial_symmetry_test(x: ndarray, y: ndarray, *, B: int = 200, seed: int = 0, alpha: float = 0.05, weights=None, bootstrap: str = 'parametric', multiplier: str = 'normal', block_length: int | str = 'auto', block_kernel: str = 'bartlett', n_jobs: int | None = None) -> RadialSymmetryResult

Test H0: the copula of (x, y) is radially symmetric (FR-10).

Parameters:

Name	Type	Description	Default
x	ndarray	pairs (raw data or pseudo-observations both work: ranks are recomputed here either way).	required
y	ndarray	pairs (raw data or pseudo-observations both work: ranks are recomputed here either way).	required
B	int	``O(B n²)`` with a per-replicate refit-free resample and rerank; ``bootstrap='multiplier'`` costs the same ``O(n²)`` once plus one ``O(n² B)`` matrix product for *all* replicates together — see the module docstring for why this is dramatically faster in practice.	200
seed	int		0
alpha	float		0.05
weights		that a caller passing ICE posteriors ``ξ`` gets an explicit error rather than a silently-wrong answer — see "Weighting" in the module docstring for why this statistic has no defensible weighted form here.	None
bootstrap	str	``'multiplier'`` (FR-10 round 2) or ``'dependent-multiplier'`` (FR-5) — see the module docstring.	'parametric'
multiplier	``'normal'`` (default) or ``'rademacher'`` — the i.i.d.	mean-0, variance-1 law of the multiplier bootstrap's ``ξ_i``. Ignored when ``bootstrap='parametric'``.	'normal'
block_length	the dependence length ``ℓ`` of the multiplier sequence,	used **only** by ``bootstrap='dependent-multiplier'``. A positive int, or ``'auto'`` (the default) for the Newey–West plug-in of :func:`pmcprg.diagnostics.dependent_multiplier.auto_block_length`. ``1`` reproduces ``bootstrap='multiplier'`` exactly.	'auto'
block_kernel	``'bartlett'`` (default) or ``'parzen'``, the smoothing	kernel of that sequence. Also only used by ``bootstrap='dependent-multiplier'``.	'bartlett'
n_jobs	int | None	replicates (audit FR-14): ``None`` (default) or ``1`` in this process, ``k > 1`` in ``k`` processes, ``-1`` one per CPU. The result is **bit-identical for every value** — a replicate draws nothing but its sampling seed, the ``b``-th draw of ``default_rng(seed)``, made here in replicate order. Starting the workers costs about a second and a replicate ``O(n²)``: at ``n = 500`` the pool is slower (0.2 s in-process, 1.0 s with 4 workers), at ``n = 2000`` it halves the time (3.6 s → 1.9 s with 4 workers, 1.8 s with 8; Apple M2 Pro, 6 performance + 4 efficiency cores, full table in the CHANGELOG entry of FR-14). Ignored by the multiplier bootstraps, which are one matrix product already.	None

Returns:

Type	Description
RadialSymmetryResult.

radial_symmetry_statistic

radial_symmetry_statistic(u: ndarray, v: ndarray) -> float

T_n of the module docstring, from pseudo-observations already in [0, 1].

Unlike :func:radial_symmetry_test, this does not rank-transform its input: pass pseudo-observations, not raw data, so that a bootstrap replicate (already uniform-margined by construction) and the observed sample are put through the identical computation.

RadialSymmetryResult dataclass

RadialSymmetryResult(statistic: float, p_value: float, reject: bool, alpha: float, n: int, tau_hat: float, B: int, n_valid: int, bootstrap: str = 'parametric', block_length: int = 1)

Typed result of :func:radial_symmetry_test.

Fields

statistic : float — the Cramér–von Mises statistic T_n (see module docstring). NaN when n < MIN_N. p_value : float — bootstrap Monte-Carlo p-value, never exactly 0 (Davison & Hinkley 1997). NaN when no replicate is available. reject : bool — p_value < alpha (H0 = "radially symmetric"). False whenever p_value is NaN. alpha : float — significance level used for reject. n : int — number of pairs the statistic was computed on. tau_hat : float — Kendall's τ of the sample. For bootstrap='parametric' it is also the calibrating Gaussian surrogate's only fitted parameter; for bootstrap='multiplier' it is reported for diagnostics only (the multiplier bootstrap fits nothing). B : int — bootstrap replicates requested. n_valid : int — replicates that produced a finite statistic. bootstrap : str — 'parametric' (default), 'multiplier' or 'dependent-multiplier', the calibration method actually used — see the module docstring. block_length : int — the multiplier dependence length ℓ actually used (FR-5). Always 1 for the two other methods, and 1 for 'dependent-multiplier' means it coincided with the i.i.d. one; report it, since the whole question FR-5 raises is how much the answer moves with ℓ.

exchangeability_test

exchangeability_test(x: ndarray, y: ndarray, *, B: int = 200, seed: int = 0, alpha: float = 0.05, weights=None, bootstrap: str = 'parametric', multiplier: str = 'normal', block_length: int | str = 'auto', block_kernel: str = 'bartlett', n_jobs: int | None = None) -> ExchangeabilityResult

Test H0: the copula of (x, y) is exchangeable, C(u,v) = C(v,u) (FR-10).

Parameters:

Name	Type	Description	Default
x	ndarray	pairs (raw data or pseudo-observations both work: ranks are recomputed here either way).	required
y	ndarray	pairs (raw data or pseudo-observations both work: ranks are recomputed here either way).	required
B	int	``O(B n²)`` with a per-replicate resample and rerank; ``bootstrap='multiplier'`` costs the same ``O(n²)`` once plus one ``O(n² B)`` matrix product for *all* replicates together — see the module docstring and :mod:`pmcprg.diagnostics.radial_symmetry`'s own docstring for why this is dramatically faster in practice.	200
seed	int		0
alpha	float		0.05
weights			None
bootstrap	str	round for this item), ``'multiplier'`` (FR-10 closing round) or ``'dependent-multiplier'`` (FR-5) — see the module docstring.	'parametric'
multiplier	``'normal'`` (default) or ``'rademacher'`` — the i.i.d.	mean-0, variance-1 law of the multiplier bootstrap's ``ξ_i``. Ignored when ``bootstrap='parametric'``.	'normal'
block_length	the dependence length ``ℓ`` of the multiplier sequence,	used **only** by ``bootstrap='dependent-multiplier'``. A positive int, or ``'auto'`` (default) for the Newey–West plug-in of :func:`pmcprg.diagnostics.dependent_multiplier.auto_block_length`. ``1`` reproduces ``bootstrap='multiplier'`` exactly.	'auto'
block_kernel	``'bartlett'`` (default) or ``'parzen'``; only used by	``bootstrap='dependent-multiplier'``.	'bartlett'
n_jobs	int | None	replicates (audit FR-14): ``None`` (default) or ``1`` in this process, ``k > 1`` in ``k`` processes, ``-1`` one per CPU. The result is **bit-identical for every value** — a replicate draws nothing but its sampling seed, the ``b``-th draw of ``default_rng(seed)``, made here in replicate order. Starting the workers costs about a second and a replicate ``O(n²)``: at ``n = 500`` the pool is slower (0.2 s in-process, 1.0 s with 4 workers), at ``n = 2000`` it halves the time (3.7 s → 1.9 s with 4 workers, 1.8 s with 8; Apple M2 Pro, 6 performance + 4 efficiency cores, full table in the CHANGELOG entry of FR-14). Ignored by the multiplier bootstraps, which are one matrix product already.	None

Returns:

Type	Description
ExchangeabilityResult.

exchangeability_statistic

exchangeability_statistic(u: ndarray, v: ndarray) -> float

T_n of the module docstring, from pseudo-observations already in [0, 1].

Unlike :func:exchangeability_test, this does not rank-transform its input: pass pseudo-observations, not raw data, so that a bootstrap replicate (already uniform-margined by construction) and the observed sample are put through the identical computation.

ExchangeabilityResult dataclass

ExchangeabilityResult(statistic: float, p_value: float, reject: bool, alpha: float, n: int, tau_hat: float, B: int, n_valid: int, bootstrap: str = 'parametric', block_length: int = 1)

Typed result of :func:exchangeability_test.

Fields

statistic : float — the Cramér–von Mises statistic T_n (see module docstring). NaN when n < MIN_N. p_value : float — parametric-bootstrap Monte-Carlo p-value, never exactly 0 (Davison & Hinkley 1997). NaN when no replicate produced a finite statistic. reject : bool — p_value < alpha (H0 = "exchangeable", i.e. C(u, v) = C(v, u)). False whenever p_value is NaN. alpha : float — significance level used for reject. n : int — number of pairs the statistic was computed on. tau_hat : float — Kendall's τ of the sample. For bootstrap='parametric' it is also the calibrating Gaussian surrogate's only fitted parameter; for bootstrap='multiplier' it is reported for diagnostics only (the multiplier bootstrap fits nothing). B : int — bootstrap replicates requested. n_valid : int — replicates that produced a finite statistic. bootstrap : str — 'parametric' (default), 'multiplier' or 'dependent-multiplier', the calibration method actually used — see the module docstring. block_length : int — the multiplier dependence length ℓ actually used (FR-5); always 1 for the two other methods.

Rosenblatt goodness-of-fit

Reduces "does (x, y) follow C_theta?" to "is (u, h(v|u)) independent-uniform on the square?", tested by a Cramér-von Mises statistic against Pi(u, v) = u*v, calibrated by a parametric bootstrap that refits theta on every replicate.

rosenblatt_gof_test

rosenblatt_gof_test(x: ndarray, y: ndarray, family_cls, *, method: str = 'tau', B: int = 200, seed: int = 0, alpha: float = 0.05, weights=None, bootstrap: str = 'parametric', multiplier: str = 'normal', block_length: int | str = 'auto', block_kernel: str = 'bartlett', n_jobs: int | None = None) -> RosenblattGoFResult

Test H0: (x, y) is drawn from family_cls (FR-10, Rosenblatt).

Parameters:

Name	Type	Description	Default
x	ndarray	pairs (raw data or pseudo-observations both work: ranks are recomputed here either way).	required
y	ndarray	pairs (raw data or pseudo-observations both work: ranks are recomputed here either way).	required
family_cls	a ``CopulaVirt`` subclass (e.g. ``CopulaGaussian``) — the	candidate family, fitted here via ``family_cls.fit``, not an already-fitted instance. ``bootstrap='parametric'`` refits it on every replicate (module docstring); ``bootstrap='multiplier'`` fits it once, on the observed sample only.	required
method	str	:meth:`CopulaVirt.fit`, used identically for the observed fit and every bootstrap refit.	'tau'
B	int	``O(B n²)``: each replicate resamples, refits (cheap: ``'tau'`` is a closed-form inversion) and recomputes the CvM statistic (``O(n²)`` for the empirical-CDF evaluation). ``bootstrap='multiplier'`` costs that ``O(n²)`` once plus one ``O(n² B)`` matrix product for *all* replicates together — see the module docstring for why this is dramatically faster in practice, and for the parameter-estimation correction it documents as omitted.	200
seed	int		0
alpha	float		0.05
weights			None
bootstrap	str	round for this item), ``'multiplier'`` (FR-10 closing round) or ``'dependent-multiplier'`` (FR-5) — see the module docstring.	'parametric'
multiplier	``'normal'`` (default) or ``'rademacher'`` — the i.i.d.	mean-0, variance-1 law of the multiplier bootstrap's ``ξ_i``. Ignored when ``bootstrap='parametric'``.	'normal'
block_length	the dependence length ``ℓ`` of the multiplier sequence,	used **only** by ``bootstrap='dependent-multiplier'``. A positive int, or ``'auto'`` (default) for the Newey–West plug-in of :func:`pmcprg.diagnostics.dependent_multiplier.auto_block_length` (selected on the *Rosenblatt-transformed* sample, which is what this statistic's empirical process is built from). ``1`` reproduces ``bootstrap='multiplier'`` exactly.	'auto'
block_kernel	``'bartlett'`` (default) or ``'parzen'``; only used by	``bootstrap='dependent-multiplier'``.	'bartlett'
n_jobs	int | None	replicates (audit FR-14): ``None`` (default) or ``1`` in this process, ``k > 1`` in ``k`` processes, ``-1`` one per CPU. The result is **bit-identical for every value** — a replicate draws nothing but its sampling seed, the ``b``-th draw of ``default_rng(seed)``, made here in replicate order. Starting the workers costs about a second; a replicate refits and transforms (``h`` once per point, in Python) on top of the ``O(n²)`` statistic. Clayton, ``n = 500``: 0.8 s in-process, 1.2 s with 4 workers (the pool is slower); ``n = 2000``: 4.2 s → 2.0 s with 4 workers, 1.9 s with 8 (Apple M2 Pro, 6 performance + 4 efficiency cores, full table in the CHANGELOG entry of FR-14). Ignored by the multiplier bootstraps, which refit nothing.	None

Returns:

Type	Description
RosenblattGoFResult.

rosenblatt_transform

rosenblatt_transform(u: ndarray, v: ndarray, copula) -> tuple[np.ndarray, np.ndarray]

Rosenblatt transform (Û, V̂) = (û, h(v̂|û)) of pseudo-observations.

Parameters:

Name	Type	Description	Default
u	ndarray	``(0, 1)`` (e.g. from :func:`_pseudo_obs` or a bootstrap replicate's own ranks — never re-derive them here so a caller controls exactly which sample the ranks come from).	required
v	ndarray	``(0, 1)`` (e.g. from :func:`_pseudo_obs` or a bootstrap replicate's own ranks — never re-derive them here so a caller controls exactly which sample the ranks come from).	required
copula	a fitted ``CopulaVirt`` instance — ``h(v|u)`` is	``copula.conditional_cdf(v, u)``.	required

Returns:

Type	Description
(U, V) : both shape ``(n,)``, clipped to ``[EPS, 1-EPS]`` (``h`` can	return exactly 0 or 1 at the boundary of a one-sided family, e.g. Clayton; the CvM statistic below only needs values inside the open square).

rosenblatt_statistic

rosenblatt_statistic(U: ndarray, V: ndarray) -> float

S_n of the module docstring, from an already-transformed sample.

S_n = n · mean_i [Ĝ_n(Û_i, V̂_i) − Û_i·V̂_i]², the Cramér–von Mises distance between the empirical CDF of (U, V) and the independent uniform reference Π(u, v) = u·v, evaluated at the sample's own points (as :func:pmcprg.diagnostics.radial_symmetry.radial_symmetry_statistic evaluates its own CvM statistic).

RosenblattGoFResult dataclass

RosenblattGoFResult(statistic: float, p_value: float, reject: bool, alpha: float, n: int, family: str, tau_hat: float, B: int, n_valid: int, bootstrap: str = 'parametric', block_length: int = 1)

Typed result of :func:rosenblatt_gof_test.

Fields

statistic : float — the Cramér–von Mises statistic S_n (see module docstring). NaN when n < MIN_N. p_value : float — parametric-bootstrap Monte-Carlo p-value, never exactly 0 (Davison & Hinkley 1997). NaN when no replicate produced a finite statistic. reject : bool — p_value < alpha (H0 = "the data follow the fitted family"). False whenever p_value is NaN. alpha : float — significance level used for reject. n : int — number of pairs the statistic was computed on. family : str — family_cls.__name__ of the candidate copula. tau_hat : float — τ of the family fitted to the observed data. B : int — bootstrap replicates requested. n_valid : int — replicates that produced a finite statistic (a fit or a Rosenblatt transform can fail on a degenerate replicate; such replicates are dropped, not counted as 0). For bootstrap='multiplier' no replicate is dropped. bootstrap : str — 'parametric' (default), 'multiplier' or 'dependent-multiplier', the calibration method actually used — see the module docstring. block_length : int — the multiplier dependence length ℓ actually used (FR-5); always 1 for the two other methods.

Pseudo-observations

Pseudo-observations and margin samples of a fitted chain, for state margins f_i and for the pair margins f_ij of a general PMC.

copula_pseudo_obs

copula_pseudo_obs(F, xi, i: int, j: int, sel=None)

Weighted pseudo-observations of the copula c_ij (DerrodePieczynski_CSDA2013 Eq. 12).

Parameters:

Name	Type	Description	Default
F			required
xi		(or one-hot draws).	required
i	int		required
j	int		required
sel	optional indices ``n`` of the transitions to keep (a thinning	stride, a window); all N−1 transitions by default.	None

Returns:

Name	Type	Description
uv	``(M, 2)`` — ``(F_ij(y_n), F_ji(y_{n+1}))``.
w	``(M,)`` — ``ξ_n(i, j)``.

margin_cdfs

margin_cdfs(model, Y, *, clip=CDF_CLIP, broadcast: bool = False) -> np.ndarray

F[n, i, j] = F_ij(y_n) for every observation and every pair.

Parameters:

Name	Type	Description	Default
model		``margin(i, j)`` (a :class:`~pmcprg.pmc.model.PMCModel`).	required
Y			required
clip		``[c, 1 − c]``, a pair ``(lo, hi)`` for ``[lo, hi]``, ``None`` for no clipping. The default is the bound the package's diagnostics used.	CDF_CLIP
broadcast	state margins only. ``False`` (default) returns a writable	``(N, K, K)`` array; ``True`` a read-only broadcast view of the ``(N, K)`` state CDFs over ``j`` — the same floats without the K copies (the GUI's form). Pair margins always return a new writable array.	False

Returns:

Type	Description
np.ndarray, shape ``(N, K, K)``, float64.
For state margins only K CDFs are evaluated and ``F[:, i, j]`` is the
vector ``F_i(Y)`` for every ``j``; for pair margins K² are evaluated.
Clipping is element-wise, so clipping the K state vectors before
broadcasting them gives the floats of clipping the broadcast array.

margin_pit_dual

margin_pit_dual(F, xi, i: int, j: int)

Weighted probability integral transform of the pair margin f_ij.

Dual view (DerrodePieczynski_CSDA2013 Eq. 12): f_ij is the left margin of the pair (i, j) and the right margin of the pair (j, i), so it accounts for

• u = F_ij(y_n) with weight ξ_n(i, j), n = 0 … N−2;
• u = F_ij(y_{n+1}) with weight ξ_n(j, i), n = 0 … N−2.

Under the model, P(x_n=i, x_{n+1}=j, y_n ∈ dy) = p_ij f_ij(y) dy and P(x_n=j, x_{n+1}=i, y_{n+1} ∈ dy) = p_ji f_ij(y) dy (DerrodePieczynski_CSDA2013 Eqs. 12–13, symmetric p), so the ξ-weighted empirical CDF of u is uniform in expectation when f_ij is right. Returned as one sample of length 2(N−1), left view first. Summed over all pairs the weights total 2(N−1): every observation is counted once as a left and once as a right observation, except the two ends.

For a state model (F[:, i, j] = F_i) the dual view of (i, j) is a subset of the state-i sample; summing its weights over j gives γ_n(i) on the left view and γ_{n+1}(i) on the right one.

Returns:

Name	Type	Description
u	``(2(N−1),)``
w	``(2(N−1),)``

margin_keys

margin_keys(model) -> list

Keys of the distinct margins of model, in a stable order.

[0, …, K−1] for state margins, [(0, 0), (0, 1), …, (K−1, K−1)] for pair margins. model.margin(*key) (or model.margin(key) for an integer key) returns the corresponding density.

margin_of

margin_of(model, key)

The margin density behind a key of :func:margin_keys.

margin_sample

margin_sample(Y, labels, key)

Observations a margin accounts for under one label path.

• key = k (state margin f_k): Y[labels == k].
• key = (i, j) (pair margin f_ij): the hard-label dual view, {y_n : (x_n, x_{n+1}) = (i, j)} followed by {y_{n+1} : (x_n, x_{n+1}) = (j, i)}. With one-hot ξ this is exactly the sample :func:margin_pit_dual weights.

labels is any path — the true labels, a posterior draw, the MPM path.

Dependent-multiplier bootstrap

Serially dependent multipliers for the bootstrap of a copula functional on a Markov chain (audit FR-5) — an i.i.d. multiplier bootstrap under-covers because consecutive pairs share an observation.

draw_multipliers

draw_multipliers(n: int, B: int, rng: Generator, law: str = 'normal', *, block_length: int = 1, kernel: str = 'bartlett') -> np.ndarray

(n, B) multipliers, mean 0, variance 1, ℓ-dependent in n.

Parameters:

Name	Type	Description	Default
n	int	``B`` columns are independent replicates, the ``n`` rows carry the serial dependence.	required
B	int	``B`` columns are independent replicates, the ``n`` rows carry the serial dependence.	required
rng	Generator	(``m = 1`` when ``block_length = 1``), so the sequence is exactly stationary rather than edge-corrected.	required
law	str	the *underlying* i.i.d. ``Z``. After smoothing with ``ℓ > 1`` the ``ξ`` are no longer ``±1`` in the Rademacher case; they remain mean 0, variance 1, which is all the multiplier CLT asks of them.	'normal'
block_length	``ℓ ≥ 1``. ``1`` returns the i.i.d. draw itself.		1
kernel	str		'bartlett'

auto_block_length

auto_block_length(u: ndarray, v: ndarray, *, kernel: str = 'bartlett', grid: tuple[float, ...] = (0.25, 0.5, 0.75), max_block: int | None = None) -> int

Newey & West (1994) plug-in ℓ for a bivariate sequence in (0,1)².

ℓ = 1 + median_{(a,b) ∈ grid²} L_NW(1{u_i ≤ a, v_i ≤ b} − C_n(a, b)) for kernel='bartlett', where L_NW is :func:pmcprg.diagnostics.model_selection.newey_west_bandwidth — the package's existing HAC bandwidth rule, reused rather than duplicated because ℓ = L + 1 is an exact identity for this kernel (module docstring). For kernel='parzen' the same target dependence length L + 1 is converted to that kernel's ℓ, whose autocorrelation reaches 2ℓ - 2: ℓ = ⌈(L + 2)/2⌉.

u, v must be in time order — the whole point is that row i is observation i of the chain. Returns at least 1, at most max_block (default n // 4).

multiplier_weights

multiplier_weights(block_length: int, kernel: str = 'bartlett') -> np.ndarray

Weights w of the moving average, normalised so Σ w_k² = 1.

block_length = 1 returns array([1.0]) for every kernel, which is what makes the dependent bootstrap reduce to the i.i.d. one exactly.

"bartlett": ℓ equal weights ℓ^{-1/2} (a rectangular window, whose self-convolution is the Bartlett kernel — module docstring). "parzen": 2ℓ-1 weights ∝ Parzen(k/ℓ), |k| ≤ ℓ-1.

multiplier_autocorrelation

multiplier_autocorrelation(block_length: int, kernel: str = 'bartlett') -> np.ndarray

[ρ(0), ρ(1), …] of the sequence, ρ(h) = Σ_k w_k w_{k+h}.

Returned for h = 0 … m-1 (m = len(w)); ρ(h) = 0 beyond. For kernel="bartlett" this is exactly 1 − h/ℓ.

MULTIPLIER_KERNELS module-attribute

MULTIPLIER_KERNELS: tuple[str, ...] = ('bartlett', 'parzen')

MULTIPLIER_LAWS module-attribute

MULTIPLIER_LAWS: tuple[str, ...] = ('normal', 'rademacher')