Skip to content

pmcprg.copulas

41 bivariate copula families (elliptical, Archimedean, survival, extreme-value, explicit, and the 90°/270° rotations for negative dependence) behind a common CopulaVirt interface (pdf, cdf, conditional_cdf, inv_h, sample, fit, plus vectorised pdf_array / logpdf_array), a bivariate-law wrapper that combines two margins with a copula, and standard-error / fit-diagnostic helpers.

Individual family classes (CopulaGaussian, CopulaClayton, …) are not listed one by one here — go through the registry (CopulaEnum) to look one up, or the "Available families" tables of the README.

Registry

CopulaEnum dataclass

CopulaEnum(ID: int, SHORT_NAME: str, LONG_NAME: str, CLASS_NAME: str, AVAILABLE: bool, PARAMETERS_SET_NAME: list[str] = list(), TAU_MIN_MAX: list[float] = list(), MODULE: str = '')

Bases: CopulaDataMixin, Enum

available classmethod

available()

Base classes

CopulaVirt

CopulaVirt(class_name: str, params: dict)

tau_range property

tau_range: tuple[float, float]

reachable_tau_bounds classmethod

reachable_tau_bounds() -> tuple[float, float] | None

(lo, hi) — the τ the family's parameter map reaches, when narrower than the registered range; None when every registered τ is reached.

A τ beyond builds the copula at the capped parameter, which then stores the τ of that parameter. Read by :meth:CopulaEnum.constructible_tau_range and :meth:CopulaEnum.reachable_tau. Default: ±reachable_tau_abs (Frank). A family whose bounds are not opposite or not known at import overrides it — Plackett, whose τ table is computed on first use and ends at −0.9935245713002141 / 0.9935245713002134.

constrain_params classmethod

constrain_params(params: dict) -> dict

Project params onto the family's admissible set.

Identity for every family whose parameters are independently bounded. Families with a joint constraint override this — BB1 needs δ < 1/(1−τ) for θ > 0, which a box cannot express. Bounded optimisers only accept per-parameter intervals, so without this hook they evaluate inadmissible points, get a constant failure penalty, and stall on that plateau with a zero gradient.

constructible_params classmethod

constructible_params(params: dict) -> dict

params itself when the constructor accepts it, else the nearest point it accepts.

params holds tau_k and any extra parameters (a TOML block's values under the constructor's names). The multistart draws move τ and each extra parameter within its own box (:meth:CopulaEnum.constructible_tau_range, EXTRA_PARAM_BOUNDS_BY_PARAM); a family whose parameters are jointly constrained overrides this hook to repair a draw that leaves the admissible set — BB1: δ < 1/(1 − τ). Unlike :meth:constrain_params, which projects onto an optimiser's fitting box, nothing is changed at a point the constructor accepts: the same dict is returned, so a draw that builds is kept bit for bit. Identity for every family with independent parameter bounds.

pdf

pdf(uv)

cdf

cdf(uv)

pdf_array

pdf_array(uv: ndarray) -> np.ndarray

Evaluate c(u, v) on M point pairs.

Parameters:

Name Type Description Default
uv (ndarray, shape(M, 2))
required

Returns:

Type Description
np.ndarray, shape (M,) — copula density values, **not floored**: 0.0
where the density underflows. ``exp(logpdf_array)`` when the subclass
provides a native log-density, the backend values otherwise.

logpdf_array

logpdf_array(uv: ndarray) -> np.ndarray

log c(u, v) on M point pairs — log(pdf_array), not floored.

The former default log(max(pdf, EPS)) capped every log-density at log EPS = −36.04 nat, flattening the likelihood surface wherever the true values are lower (−289 to −3086 nat for Student at τ = 0.95, ν = 30; audit RB-4). An exact zero now gives −∞, which the likelihood objectives handle (:func:pmcprg.copulas._fit._weighted_log_density_sum).

cdf_array

cdf_array(uv: ndarray) -> np.ndarray

Evaluate C(u, v) on M point pairs.

Same fast-path strategy as :meth:pdf_array: use the optional vectorised backend self._model if present (its .cdf accepts a (M, 2) ndarray and returns an (M,) ndarray), otherwise fall back to a Python loop.

Subclasses with closed-form CDFs (FGM, AMH, Joe, …) can override this with a vectorised expression for ~50× speed-ups on bootstrap GoF.

Parameters:

Name Type Description Default
uv (ndarray, shape(M, 2))
required

Returns:

Type Description
np.ndarray, shape (M,) — clipped to [0, 1]; NaN where the backend
returns a non-finite value.
A non-finite backend value used to be replaced by EPS, a plausible
number in the lower corner and a wrong one wherever C is close to 1
(audit RB-9). NaN makes the failure visible to the caller instead.

fit classmethod

fit(data: ndarray, method: str = 'tau', *, weights=None, pseudo_obs: bool = False) -> FitResult

Fit a copula's free parameters from a (n, 2) data array.

Raw data is rank-transformed to pseudo-observations û = rank/(n+1) before estimation — no distributional assumption on the margins.

Parameters:

Name Type Description Default
data ndarray
required
method (tau, mle)

'tau' — inversion of Kendall's τ (itau; Genest & Rivest 1993), O(n log n), the default. A family with more parameters than τ gets the others by maximum likelihood with τ held at τ̂ — a profile likelihood, VineCopula's rule for Student's ν (FR-12; before, these families silently ran the joint MLE and reported method='mle'). 'mle' — maximise ∑ w log c(û_i, v̂_i). For 1-parameter families a 1-D Brent scalar search on the padded τ-range; for two-parameter families that do not override fit the joint optimiser of the ICE M-step (:func:pmcprg.copulas._fit._fit_two_parameter_mle). On rank pseudo-observations this is the pseudo-maximum-likelihood estimator of Genest, Ghoudi & Rivest (1995).

'tau'
weights (array - like, shape(n))

Observation weights (FR-12) — frequency weights, the package's convention: an integer weight counts its row that many times. Validated by :func:pmcprg.copulas._fit.validate_weights (1-D, one per row, finite, ≥ 0, positive sum; ValueError otherwise). Rows of weight zero are dropped; if every other weight is 1 the result is the unweighted fit of those rows, bit for bit (so {0, 1} weights give the fit on the subset). Otherwise the weighted engine of ICE's M-step (:mod:pmcprg.copulas._weighted) is used:

  • pseudo-observations by the weighted empirical CDF Σ_k w_k 1{x_k ≤ x}/(Σw + 1) (ICE's copula_margins = "empirical" convention, audit FR-7 a) unless pseudo_obs;
  • 'tau' — the weighted Kendall τ, a weighted τ-b as in pyvinecopulib and VineCopula (the τ-a when no two points of positive weight tie), then as above with the weighted profile;
  • 'mle' — the weighted MLE of ICE's M-step: a Brent search on the Σw-normalised objective (xatol 10⁻⁶) for one parameter, the joint optimiser with gtol scaled by the mean weight for more — parameters invariant to the weights' scale.

The log-likelihood is the total Σ wᵢ log cᵢ and the BIC charges k·log Σw (see :class:FitResult).

None
pseudo_obs bool

True when data already are pseudo-observations in (0, 1)²: they are used as they are (ValueError if a value is outside the open square or not finite). False rank-transforms them.

False

Returns:

Type Description
FitResult
Notes

With 'tau', τ̂ is clipped into the constructible τ-range: a singular registered endpoint (|τ| = 1) is replaced by the padded value next to it, an admissible endpoint is kept (τ̂ ≤ 0 → ε for Clayton/GH/Joe), and a τ̂ beyond the family's reachable |τ| goes to that bound (Frank). With 'mle', a τ̂ found beyond the reachable |τ| — where the likelihood is flat — is returned as the bound too, so tau_k is always the τ the fitted copula uses. FitResult.converged (and failed, message) report whether the optimiser converged; FitResult.diagnostics the gradient, Hessian eigenvalues and boundary flag of the estimate.

References
  • Genest, C. & Rivest, L.-P. (1993). Statistical inference procedures for bivariate Archimedean copulas. JASA 88(423), 1034–1043.
  • Genest, C., Ghoudi, K. & Rivest, L.-P. (1995). A semiparametric estimation procedure of dependence parameters in multivariate families of distributions. Biometrika 82(3), 543–552.

fit_best staticmethod

fit_best(data: ndarray, families: list | None = None, method: str = 'tau', *, weights=None, criterion: str = 'aic', pseudo_obs: bool = False) -> FitBestResults

Fit each candidate family and rank the comparable fits by criterion.

Parameters:

Name Type Description Default
data ndarray
required
families list | None
     family except Product.
None
method str
'tau'
weights
     (see :meth:`fit`; FR-12).
None
criterion str
     ranking, unchanged), ``'bic'`` (ascending BIC: ``k·log n``
     unweighted, ``k·log Σw`` weighted, as ICE's ``bic``
     selection criterion — VineCopula's ``BiCopSelect`` uses
     ``log n`` rows) or ``'loglik'`` (descending
     log-likelihood, ICE's ``mle`` criterion). With weights,
     ``method='mle'`` and a criterion, the ranking is ICE's
     weighted family selection
     (:func:`pmcprg.copulas._weighted._select_and_fit_copula`)
     over the same candidates.
'aic'
pseudo_obs passed to each ``cls.fit`` (see :meth:`fit`).
False

Returns:

Type Description
FitBestResults

A list[FitResult] sorted by the criterion (best first), holding only the fits obtained with the requested method.

Notes

Two attributes of the result report what the ranking leaves out (audit RB-10):

  • other_method — fits a family could only produce with another method (not ranked: their log-likelihood is not comparable). Since FR-12 every registered family fits by either method — Student, BB1 and the other multi-parameter families answer 'tau' by itau with a profile MLE of their other parameters — so this is empty unless a family overrides fit.
  • failures — (class name, "ExceptionType: message") for every family whose fit raised (formerly logged and dropped).

Both are also logged at WARNING. converged on each result says whether its likelihood optimiser converged.

conditional_cdf

conditional_cdf(v: float, u: float) -> float

h(v|u) = ∂C(u,v)/∂u via central finite differences. Subclasses override analytically.

transposed

transposed() -> CopulaVirt

The copula of (V, U) when (U, V) has this one: Cᵀ(u, v) = C(v, u).

The h-functions condition on the first argument; conditioning on the second is conditioning on the first of the transpose:

∂C/∂v (u, v) = P(U ≤ u | V = v) = ``transposed().conditional_cdf(u, v)``,

inverted in u by transposed().inv_h(w, v). self for an exchangeable family. Checked for every registered family (pmcprg/tests/test_copula_transpose.py), and against the h2 and h⁻¹ of pyvinecopulib, VineCopula and R copula (pmcprg/tests/test_parity_interior.py).

inv_h

inv_h(w: float, u: float) -> float

Return v such that h(v|u) = w (Rosenblatt inverse step).

Default numerical implementation: bracket h(v|u) - w on (EPS, 1-EPS) and solve via Brent's method. Brent requires the function to be monotone in v on that bracket — this holds for any valid copula (h is a conditional CDF in v) but the assertion below guards against an ill-defined custom copula.

Subclasses with closed-form inverses (Gaussian, Clayton, Frank, …) override this for ~50× speed-up.

Saturation (kept deliberately, audit RB-9): when w lies outside the range [h(EPS|u), h(1−EPS|u)] that the bracket can attain, the root is outside (EPS, 1−EPS) and the method returns the nearer bracket end, EPS or 1−EPS. This is the value of the exact inverse clamped to the grid every sampler in the package works on — not a failure — so it is not reported.

inv_h_array

inv_h_array(w: ndarray, u: ndarray) -> np.ndarray

Vectorised version of :meth:inv_h.

Default fallback: scalar Python loop. Subclasses with closed-form inv_h (Gaussian, Clayton, Frank) override this with a vectorised expression for ~10× speed-up on bulk sampling.

Parameters:

Name Type Description Default
w np.ndarray, shape (M,) — uniform variates.
required
u np.ndarray, shape (M,) — conditioning values.
required

Returns:

Type Description
np.ndarray, shape (M,) — v values such that h(v_i | u_i) = w_i.

tail_dependence

tail_dependence() -> tuple[float, float]

Lower/upper tail dependence (λ_L, λ_U).

Default numerical implementation via diagonal CDF limits: λ_L = lim_{u→0+} C(u,u)/u λ_U = lim_{u→1−} (1 − 2u + C(u,u))/(1 − u) Subclasses with analytical formulas override this. Returns (nan, nan) if the copula has no closed-form CDF.

standard_errors

standard_errors(uv: ndarray, weights: ndarray | None = None, method: str = 'mle', *, ranks: bool = True)

Asymptotic standard errors of this copula's parameters as estimated on uv.

Opt-in and side-effect free: nothing in :meth:fit changes. The copula's current parameters are taken as the estimate; uv are the pseudo-observations it was computed from, weights optional frequency weights (n_eff = Σw).

  • method='mle' — sandwich of the pseudo-likelihood with the estimated-rank corrections (Genest, Ghoudi & Rivest 1995, doi:10.1093/biomet/82.3.543); ranks=False for known margins.
  • method='tau' — 16·Var{2C(U,V) − U − V}/n for Kendall's τ̂, carried to θ by dθ/dτ (Genest & Favre 2007, doi:10.1061/(ASCE)1084-0699(2007)12:4(347); Kojadinovic & Yan 2010, doi:10.1016/j.insmatheco.2010.03.008). One-parameter families only.

Returns a :class:pmcprg.copulas._stderr.StandardErrors (estimate, se, cov, ci(level), at_boundary, n_eff). At a boundary of the parameter space (independence end of Clayton/GH/Joe, BB1 δ = 1, Student ν at its bound) at_boundary is set and no Wald interval is reported (Self & Liang 1987, doi:10.1080/01621459.1987.10478472). See :mod:pmcprg.copulas._stderr for the formulas and their limits (i.i.d. pairs; FR-4 inside ICE and FR-5 not covered).

sample

sample(n: int = 500, seed: int | None = None) -> np.ndarray

Draw n samples from the copula on [0,1]² via Rosenblatt inversion.

For each sample::

u ~ Uniform(EPS, 1-EPS)
w ~ Uniform(EPS, 1-EPS)
v = inv_h(w, u)

Subclasses with a closed-form inv_h_array (Gaussian, Clayton, Frank) get a single vectorised call for all n samples; others fall back to per-sample Brent inversion via inv_h.

FitResult dataclass

FitResult(copula: 'CopulaVirt', method: str, tau_k: float, log_likelihood: float, n_obs: int, uv: ndarray, converged: bool = True, message: str = '', n_iter: int | None = None, n_eval: int | None = None, weights: ndarray | None = None)

Returned by :meth:CopulaVirt.fit.

Fields

copula, method, tau_k : the fitted copula, the method ('tau' or 'mle') and its τ. log_likelihood : Σᵢ wᵢ log c(ûᵢ, v̂ᵢ) at the estimate (wᵢ = 1 unweighted) — the total, not divided by Σw, as in VineCopula's weighted MLE and in ICE's selection scores. pyvinecopulib reports (n/Σw)·Σ wᵢ log cᵢ instead (its weights rescaled to mean 1); VineCopula's itau ignores the weights. n_obs : number of observations of positive weight. uv : the pseudo-observations the estimate was computed from, shape (n_obs, 2). converged : the optimiser's verdict: False when the search behind the estimate did not converge (the copula is then the best point it found, and a WARNING was logged). True for a moment estimate (a one-parameter 'tau' fit). message : what the optimiser said ("" when not recorded). n_iter, n_eval : its iterations and log-likelihood evaluations (0 for a moment estimate, None when not recorded). weights : the weights of the n_obs rows for a weighted fit, None for an unweighted one (unit weights are fitted as unweighted, :meth:CopulaVirt.fit).

failed (not converged) is the public counterpart of ICE's internal FIT_FAILED_KEY, and :attr:diagnostics the per-fit convergence diagnostics (gradient, Hessian eigenvalues, boundary flag, counts — :mod:pmcprg.copulas._fit_diagnostics), computed on first access only.

Information criteria: k free parameters and n_eff = :attr:n_eff (n_obs unweighted, Σw weighted, the size ICE's BIC charges): AIC 2k − 2ℓ, BIC k·log n_eff − 2ℓ (floored at n_eff = 1 when weighted, as in ICE), AICc and HQC with n_eff for n. VineCopula's BIC charges log n (rows) instead of log Σw.

failed property

failed: bool

not converged — the estimate is the best point found, not a certified optimum.

weighted property

weighted: bool

True for a weighted fit (weights given, not all ones).

n_eff property

n_eff: float

Effective sample size: n_obs unweighted, Σw weighted (frequency weights).

aic property

aic: float

Akaike Information Criterion: 2k − 2·loglik (smaller is better).

bic property

bic: float

Bayesian Information Criterion: k·log(n) − 2·loglik (smaller is better).

n is n_obs unweighted and max(Σw, 1) weighted — the penalty of ICE's bic selection score (pmcprg.pmc.ice), so a weighted fit_best(criterion='bic') ranks as ICE selects.

aicc property

aicc: float

Corrected AIC for small samples: AIC + 2k(k+1)/(n−k−1), n = n_eff.

hqc property

hqc: float

Hannan-Quinn: 2k·log(log n) − 2·loglik, n = n_eff (asymptotically less biased than BIC).

diagnostics cached property

diagnostics

Convergence diagnostics of this fit — :class:pmcprg.copulas.FitDiagnostics.

The gradient of the (weighted) log-likelihood at the estimate on the parameter scale, the eigenvalues of its Hessian (negative at an interior maximum), a boundary flag and the optimiser's counts, as GJRM's conv.check() reports them. Computed on first access (a few likelihood evaluations: 3 for one parameter, 9 for two), then cached; nothing is computed during the fit. See :mod:pmcprg.copulas._fit_diagnostics.

gof_test

gof_test(B: int = 100, seed: int | None = None, *, n_jobs: int | None = None) -> GoFResult

Cramér-von Mises GoF test via parametric bootstrap.

H₀: data was generated by the fitted copula family. Unweighted fits only (NotImplementedError for a weighted one). Returns S_n and a bootstrap p-value over B replicates. The copula is re-fitted on every replicate, so this is the parametric bootstrap of Genest & Rémillard (2008, Ann. IHP 44(6):1096–1127) for the composite null, as reviewed by Genest, Rémillard & Beaudoin (2009, Insurance Math. Econom. 44(2):199–213). The p-value is (1 + #{S*_b ≥ S_n}) / (n_valid + 1) (Davison & Hinkley 1997, ch. 4) — never 0; GRB 2009 (App. A) write the plain proportion, which differs by at most 1/(B + 1).

n_jobs (audit FR-14) — None (default) or 1 runs the replicates in this process, k > 1 in k worker processes, -1 in 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), and those draws are made here, in replicate order, whatever n_jobs is. Starting the workers costs about a second, and a replicate is one sample and one refit — a millisecond or two for most families, so the pool is slower (Clayton by τ, n = 500, B = 100: 0.12 s in-process, 1.0 s with 4 workers); tens of ms when sample inverts h numerically, where it pays (BB1 by MLE: 4.5 s → 2.2 s with 4 workers, 2.0 s with 8). Measured on an Apple M2 Pro (6 performance + 4 efficiency cores); the full table is in the CHANGELOG entry of FR-14. Warnings and log records raised in the workers are re-emitted here (:mod:pmcprg._parallel).

bootstrap_ci

bootstrap_ci(B: int = 500, alpha: float = 0.05, seed: int | None = None, *, n_jobs: int | None = None) -> tuple[float, float]

Non-parametric bootstrap percentile CI for tau_k at level (1-alpha).

Resamples the pseudo-observations with replacement B times, refits with the same method, and returns the (alpha/2, 1-alpha/2) quantiles. Unweighted fits only (NotImplementedError for a weighted one).

n_jobs (audit FR-14) — None (default) or 1 runs the replicates in this process, k > 1 in k worker processes, -1 in one per CPU. The interval is bit-identical for every value: the resampling indices of replicate b are the b-th draw of default_rng(seed), made here, in replicate order, and sent to the worker with the replicate (at most 2 · n_jobs chunks of them in flight). Starting the workers costs about a second, and a replicate is one refit — under a millisecond for a moment fit, so the pool is slower (Clayton by τ, n = 500, B = 500: 0.25 s in-process, 1.0 s with 4 workers); tens of ms for a joint MLE of several parameters, where it pays (Student by MLE: 22.7 s → 6.9 s with 4 workers, 4.9 s with 8). Measured on an Apple M2 Pro (6 performance + 4 efficiency cores); the full table is in the CHANGELOG entry of FR-14. Warnings and log records raised in the workers are re-emitted here (:mod:pmcprg._parallel).

standard_errors

standard_errors(*, ranks: bool = True)

Asymptotic standard errors of this fit — opt-in, nothing is cached.

copula.standard_errors(self.uv, method=self.method, ranks=ranks): the pseudo-likelihood sandwich with the estimated-rank corrections for method='mle' (Genest, Ghoudi & Rivest 1995), the variance of the inversion of Kendall's τ for method='tau' (Genest & Favre 2007; Kojadinovic & Yan 2010). Returns a :class:pmcprg.copulas._stderr.StandardErrors; see :mod:pmcprg.copulas._stderr. Unlike :meth:bootstrap_ci, it is analytic (no refit) and reports SEs on τ and on the native parameters. A weighted fit passes its weights (frequency weights, n_eff = Σw).

cv_loglik

cv_loglik(K: int = 5, seed: int | None = None) -> float

K-fold CV log-likelihood (held-out test loglik summed over folds).

For each fold, the copula is refitted on the K−1 training folds (using the same method) and its log-density is summed over the held-out test fold's pseudo-observations. Higher is better. Unweighted fits only (NotImplementedError for a weighted one).

plot_diagnostics

plot_diagnostics(plot_dir: str, prefix: str = '') -> None

6-panel diagnostic plot: (0,0) pseudo-observations + fitted PDF contours (0,1) PP plot C_n vs C_θ at the data points (0,2) lower tail dependence λ̂_L(u) vs fitted λ_L (1,0) empirical copula heatmap C_n(u,v) (1,1) residuals heatmap C_n − C_θ (1,2) upper tail dependence λ̂_U(u) vs fitted λ_U

Unweighted fits only: the empirical copula and tail curves are the unweighted ones (NotImplementedError for a weighted fit).

FitDiagnostics dataclass

FitDiagnostics(names: tuple, estimate: ndarray, gradient: ndarray, hessian: ndarray, eigenvalues: ndarray, at_boundary: bool, boundary: tuple, converged: bool, message: str, n_iter: int | None, n_eval: int | None, steps: ndarray)

Convergence diagnostics of one fit — see the module docstring.

Fields

names : the free parameters, on the constructor's scale (('tau_k',), ('tau_k', 'df'), …); empty for the Product copula. estimate : their fitted values. gradient : ∂ℓ/∂p at the estimate (ℓ the weighted log-likelihood). hessian : ∂²ℓ/∂p∂pᵀ at the estimate (NaN entries where no difference could be taken). eigenvalues : eigenvalues of hessian, ascending (NaN if it has a NaN). at_boundary : the estimate is on, or within the documented tolerance of, the boundary of its admissible range (module docstring). boundary : one sentence per boundary condition met. converged : the optimiser's own verdict (FitResult.converged). message : the optimiser's message (FitResult.message). n_iter : optimiser iterations (None when not reported; 0 for a moment estimate). n_eval : log-likelihood evaluations by the optimiser (same). steps : the finite-difference step of each coordinate.

negative_definite property

negative_definite: bool

Every eigenvalue of the Hessian is finite and < 0 — an interior maximum.

max_abs_gradient property

max_abs_gradient: float

max |∂ℓ/∂p| — GJRM's conv.check() headline number (NaN if undefined).

newton_decrement property

newton_decrement: float

½ gᵀ(−H)⁻¹g in nats: the log-likelihood a Newton step would still gain.

Invariant under reparametrisation. NaN unless the Hessian is negative definite.

summary

summary() -> str

A few lines in the manner of GJRM's conv.check().

Weighted fitting

fit and fit_best accept observation weights (weights=); ICE's M-step uses the same engine (pmcprg.copulas._weighted).

validate_weights

validate_weights(weights, n: int | None = None, *, allow_zero_sum: bool = False) -> np.ndarray

weights as a float array, after the package's one check of observation weights.

Every weighted computation of the package validates its weights here — the public CopulaVirt.fit(weights=…) and fit_best, the copula step of ICE/SEM, the weighted log-likelihood behind every objective (:func:_weighted_log_density_sum), the standard errors, the robust fits and the model-selection tests — so they accept and refuse the same weights with the same messages. A valid weight vector is

  • one-dimensional,
  • of length n (one weight per observation) when n is given,
  • finite and non-negative (a zero weight drops its observation),
  • of positive sum — unless allow_zero_sum, for the low-level sums whose callers handle an all-zero vector themselves (ICE skips a pair of states with Σw < 1e-12 before it gets here).

Raises ValueError naming the first condition that fails. Returns np.asarray(weights, dtype=float) — the same array, not a copy, when it already is a float array — so validating costs no copy and changes no value.

weighted_kendall_tau

weighted_kendall_tau(u, v, weights=None) -> float

Weighted Kendall's τ of two samples — the coefficient the weighted fit inverts.

The weighted τ-b of pyvinecopulib and VineCopula, which is the weighted τ-a when no two points of positive weight tie (the engine's _weighted_kendall_tau, O(n log n)); weights=None gives the τ-b of scipy.stats.kendalltau. Only ranks matter: u, v may be raw data or pseudo-observations. Weights are validated (:func:validate_weights); u and v must have one value per weight.

weighted_pseudo_obs

weighted_pseudo_obs(data, weights) -> np.ndarray

Weighted pseudo-observations of an (n, 2) sample: each column through its weighted ECDF.

û_ij = Σ_k w_k 1{x_kj ≤ x_ij} / (Σ_k w_k + 1) — the convention of ICE's posterior-weighted empirical margins (audit FR-7 a, :func:_weighted_ecdf) — clipped to [EPS, 1 − EPS] as those are.

  • Unit weights give rank/(n + 1), the rank pseudo-observations of the unweighted fit, exactly when there are no ties (ties share the upper end of their block here, the mid-rank in the unweighted fit).
  • {0, 1} weights give the rank pseudo-observations of the kept subset, since a zero weight adds nothing to any F̂.
  • The + 1 is on the scale of the weights: they are frequency weights (an integer weight counts its observation that many times, the package's convention, as for n_eff = Σw and the BIC), so multiplying every weight by c ≠ 1 moves the pseudo-observations by O(1/Σw) (measured: 4.7·10⁻³ for ×10³ at Σw = 210). Pass pseudo-observations with pseudo_obs=True for a fit that does not depend on the weights' scale. Weights summing to less than 10 put every pseudo-observation below Σw/(Σw + 1) < 0.91 and log a WARNING — probability weights, summing to 1, put them all below ½.

Raises ValueError for data that are not (n, 2) and finite, and for invalid weights (:func:validate_weights).

GoFResult dataclass

GoFResult(statistic: float, p_value: float, B: int, n_valid_bootstrap: int, bootstrap_stats: ndarray)

Returned by :meth:FitResult.gof_test — Cramér-von Mises GoF test.

Bivariate joint laws

Combine two margins (state or pair densities) with a copula into the joint law f_ij(y_1, y_2) = f_ij(y_1) f_ji(y_2) c_ij(F_ij(y_1), F_ji(y_2)) (DerrodePieczynski_CSDA2013 Eq. 12).

BivariateLaw

BivariateLaw(copula, left_margin: tuple, right_margin: tuple, quantile_range: tuple[float, float] = (0.05, 0.95))

Bivariate distribution built from a copula and two marginals (Sklar's theorem):

f(x, y) = f1(x) · f2(y) · c(F1(x), F2(y))

Parameters:

Name Type Description Default
copula
required
left_margin tuple
required
right_margin tuple
required
quantile_range tuple[float, float]
(0.05, 0.95)

log_pdf

log_pdf(xy) -> float

log f(x,y) = log f1(x) + log f2(y) + log c(F1(x), F2(y)).

pdf

pdf(xy) -> float

Joint PDF at (x, y) in the original scale.

cdf

cdf(xy) -> float

Joint CDF: C(F1(x), F2(y)).

conditional_log_pdf

conditional_log_pdf(y_cond: float, y_obs: float, which: str = 'left') -> float

log p(y_cond | y_obs) = log f_cond(y_cond) + log c(u_left, u_right).

which='left' → condition on left margin value, evaluate for right which='right' → condition on right margin value, evaluate for left

conditional_pdf

conditional_pdf(y_cond: float, y_obs: float, which: str = 'left') -> float

p(y_cond | y_obs) = f_cond(y_cond) · c(u_left, u_right).

conditional_cdf

conditional_cdf(y_cond: float, y_obs: float, which: str = 'left') -> float

P(Y_cond ≤ y_cond | Y_obs = y_obs).

∂C/∂u at (u_obs, u_cond) — the copula's h-function conditional_cdf — when the left margin is observed; ∂C/∂v at (u_cond, u_obs), the h-function of the transposed copula, when the right one is.

conditional_law

conditional_law(y_obs: float, which: str = 'left') -> ConditionalLaw

Return a frozen conditional distribution p(Y_cond | Y_obs = y_obs).

which='left' → condition on left margin, sample/evaluate the right which='right' → condition on right margin, sample/evaluate the left

sample

sample(n: int = 1) -> np.ndarray

Draw n iid samples from the joint distribution. Returns shape (n, 2).

Vectorised Rosenblatt: y_left ~ F_left (n at once via rvs(size=n)) u_left = F_left(y_left) w ~ Uniform(0, 1) (n at once) v = inv_h_array(w, u_left) y_right = F_right⁻¹(v)

sample_conditional

sample_conditional(y_obs: float, which: str = 'left', n: int = 1) -> np.ndarray

Draw n iid samples from p(Y_cond | Y_obs = y_obs). Returns shape (n,).

Vectorised Rosenblatt — scales as O(n) numpy ops, not O(n) Python loop iterations.

fit classmethod

fit(data: ndarray, copula_class, left_family, right_family, copula_method: str = 'tau', quantile_range: tuple[float, float] = (0.05, 0.95)) -> BivariateFitResult

Two-step IFM fit: marginal MLE, then copula fit on the data.

Step 1: each marginal distribution is MLE-fitted to its column. Step 2: the copula is fitted on the data (rank-based pseudo-obs, see :meth:CopulaVirt.fit).

Parameters:

Name Type Description Default
data ndarray
required
copula_class
required
left_family
required
right_family
required
copula_method str
'tau'
quantile_range forwarded to BivariateLaw plotting grid
(0.05, 0.95)

Returns:

Type Description
BivariateFitResult

fit_best staticmethod

fit_best(data: ndarray, left_family, right_family, copula_families: list | None = None, copula_method: str = 'tau', quantile_range: tuple[float, float] = (0.05, 0.95)) -> list[BivariateFitResult]

Fit each candidate copula family with given margins, sorted by AIC.

Parameters:

Name Type Description Default
data ndarray
required
left_family
required
right_family
required
copula_families list of CopulaVirt subclasses, or None for a default set
          (all available 1-parameter copulas except Product)
None
copula_method str
'tau'

Returns:

Type Description
list[BivariateFitResult] sorted by AIC ascending. Failures are logged.

ConditionalLaw

ConditionalLaw(parent: BivariateLaw, y_obs: float, which: str)

Frozen conditional distribution p(Y_cond | Y_obs = y_obs).

Obtained via BivariateLaw.conditional_law(y_obs, which). All methods delegate to the parent BivariateLaw and share its RNG.

Parameters:

Name Type Description Default
parent BivariateLaw
required
y_obs float
required
which str
 'right' → condition on right margin, output is left
required

pdf

pdf(y: float) -> float

Conditional density p(y | y_obs).

log_pdf

log_pdf(y: float) -> float

Log conditional density.

cdf

cdf(y: float) -> float

Conditional CDF P(Y_cond ≤ y | y_obs).

sample

sample(n: int = 1) -> np.ndarray

Draw n iid samples. Returns shape (n,).

BivariateFitResult dataclass

BivariateFitResult(bivariate: 'BivariateLaw', copula_fit: FitResult, left_params: tuple, right_params: tuple, log_likelihood: float, n_obs: int, data: ndarray)

Returned by :meth:BivariateLaw.fit.

aicc property

aicc: float

Corrected AIC for small samples: AIC + 2k(k+1)/(n−k−1).

hqc property

hqc: float

Hannan-Quinn: 2k·log(log n) − 2·loglik.

gof_test

gof_test(B: int = 100, seed: int | None = None) -> GoFResult

Cramér-von Mises GoF test on the joint distribution.

Tests H₀: data was generated by the fitted bivariate (margins + copula). Bootstrap resamples from the fitted BivariateLaw and refits at each iteration — heavier than the copula-only test in :class:FitResult, but penalises misspecified margins as well.

bootstrap_ci

bootstrap_ci(B: int = 500, alpha: float = 0.05, seed: int | None = None) -> BivariateBootstrapCI

Non-parametric bootstrap percentile CI on every fitted parameter.

Resamples (X_i, Y_i) pairs with replacement, refits the entire BivariateLaw (margins + copula), and returns (alpha/2, 1-alpha/2) quantiles for τ_k and each marginal parameter.

cv_loglik

cv_loglik(K: int = 5, seed: int | None = None) -> float

K-fold CV log-likelihood for the full bivariate model.

For each fold, both margins (MLE) and the copula are refitted on the training set, and ∑ log f(x,y) is accumulated on the held-out test set. Higher is better. Cleaner than AIC: penalises overfitting honestly.

Numerical floor

When a fitted margin's support excludes some test point (typical for expon / triang / similar location-shifted families on small-N folds), log f(x, y) = -∞. We clip these contributions at -100 nats rather than log(MIN_POSITIVE) ≈ -708. The choice is intentional: * -100 keeps the CV score finite and on a comparable scale across folds (a single -708 outlier would dominate the sum); * -100 still penalises support violations heavily (e^{-100} ≈ 0, a fold with several violations clearly loses to its peers); * the absolute value of the CV score is rarely used — only the relative ranking across competing models matters.

plot_diagnostics

plot_diagnostics(plot_dir: str, prefix: str = '') -> None

6-panel diagnostic plot: (0,0) data scatter + fitted joint PDF contours (0,1) QQ plot — left margin (0,2) lower tail dependence λ̂_L(u) vs fitted λ_L (on parametric pseudo-obs) (1,0) QQ plot — right margin (1,1) PP plot of the joint CDF F_n vs F_θ (1,2) upper tail dependence λ̂_U(u) vs fitted λ_U

BivariateBootstrapCI dataclass

BivariateBootstrapCI(tau_k: tuple[float, float], left_params: list[tuple[float, float]], right_params: list[tuple[float, float]], alpha: float, B: int, n_valid: int)

Returned by :meth:BivariateFitResult.bootstrap_ci.

Fitting diagnostics

Standard errors, and the independence / sub-model likelihood-ratio tests (audit FR-4).

standard_errors

standard_errors(copula, uv, weights=None, method: str = 'mle', *, ranks: bool = True) -> StandardErrors

Asymptotic standard errors of the parameters of copula estimated on uv.

Parameters:

Name Type Description Default
copula
required
uv
  computed from** (``FitResult.uv`` after :meth:`CopulaVirt.fit`).
required
weights optional (n,) non-negative frequency weights (module docstring).
None
method str
  1995); ``'tau'`` — variance of the inversion of Kendall's τ
  (Genest & Favre 2007; Kojadinovic & Yan 2010). ``'tau'`` needs a
  one-parameter family: τ does not identify Student's ν or BB1's δ.
'mle'
ranks bool
  pseudo-observations, whose estimated margins add the W terms;
  ``False`` for margins taken as known.
True

Returns:

Type Description
StandardErrors

Raises:

Type Description
ValueError for a family without a free parameter, an unknown

method, 'tau' on a two-parameter family, or invalid uv / weights.

NotImplementedError for a multi-parameter family other than Student and BB1.

StandardErrors dataclass

StandardErrors(family: str, method: str, ranks: bool, names: tuple, estimate: dict, se: dict, cov: ndarray, at_boundary: bool, boundary: tuple, n_obs: int, n_eff: float)

Asymptotic standard errors of a copula estimate — see the module docstring.

Fields

family : class name of the copula. method : 'mle' (pseudo-likelihood sandwich) or 'tau' (τ inversion). ranks : whether the estimated-margin correction (W terms) was applied (always True for 'tau', whose variance already is the rank statistic's). names : reported quantities, in the order of cov: ('tau_k', 'theta') for one-parameter families (theta is the family's native parameter — ρ for the Gaussian), ('tau_k', 'df', 'rho') for Student, ('tau_k', 'delta', 'theta') for BB1. estimate : {name: value} at which the variance was evaluated. se : {name: standard error}; NaN where it could not be computed, and for 'mle' at a boundary (module docstring). cov : covariance matrix of names (rank = number of free parameters; the native parameters are deterministic functions of the free ones). at_boundary : the estimate sits on the boundary of the parameter space (see boundary); Wald intervals are not reported. boundary : one sentence per boundary condition met (empty if interior). n_obs : number of positive-weight observations. n_eff : Σ w — the sample size the variance refers to.

ci

ci(level: float = 0.95, name: str = 'tau_k') -> tuple[float, float]

Wald interval estimate ± z_{(1+level)/2} · se for name.

(nan, nan) when the estimate is at a boundary of the parameter space — the estimator is then not asymptotically normal (Self & Liang 1987) — or when the standard error is not finite. The interval is on the natural scale and is not clipped to the parameter space.

cov_of

cov_of(*names: str) -> np.ndarray

Covariance sub-matrix of the requested quantities, in the given order.

independence_lr_test

independence_lr_test(family, uv, weights=None) -> IndependenceLRTest

Likelihood-ratio test of independence within a one-parameter family.

The (weighted) pseudo-likelihood ℓ(τ) = Σ w log c(û, v̂; τ) is maximised over the padded τ-range, as :meth:CopulaVirt.fit does, and compared with the independence copula, ℓ(Π) = 0: LR = 2 max(ℓ(τ̂), 0).

  • Independence inside the family (Gaussian, Frank, FGM, AMH, Plackett): LR → χ²₁.
  • Independence at an end of the parameter space (Clayton, Gumbel, Joe, their survival versions, Cubic Section — τ ≥ 0 only — and the 90°/270° rotations of FR-8 whose range ends at τ = −ε, e.g. CopulaClayton90 — τ ≤ 0 only): the restricted maximiser sits on the boundary about half of the time under H₀, and LR → ½χ²₀ + ½χ²₁ (Self & Liang 1987): p = ½ P(χ²₁ ≥ LR) for LR > 0, p = 1 for LR = 0.

On rank-based pseudo-observations the W corrections of the pseudo- likelihood score vanish under independence (∫ φ(u, v) dv = 0 for every u when c ≡ 1), so these reference distributions hold unchanged (Genest, Ghoudi & Rivest 1995). Weights are frequency weights; with fractional weights ≤ 1 the statistic is stochastically smaller than its reference and the test conservative.

Parameters:

Name Type Description Default
family
required
uv
required
weights optional (n,) non-negative weights.
None

Raises:

Type Description
ValueError when the family does not contain independence (A12, A14 and
their rotations, whose ranges stop at ±1/3), has
no free parameter (Product) or has two parameters.

IndependenceLRTest dataclass

IndependenceLRTest(family: str, statistic: float, p_value: float, null_distribution: str, boundary: bool, tau_k: float, log_likelihood: float, n_obs: int, n_eff: float)

Likelihood-ratio test of independence in a one-parameter family.

Fields

family : class name. statistic : 2 (ℓ(τ̂) − ℓ(Π)) with ℓ = Σ w log c and ℓ(Π) = 0; ≥ 0. p_value : from null_distribution. null_distribution : 'chi2(1)' when independence (τ = 0) is interior to the family, '0.5*chi2(0) + 0.5*chi2(1)' when it is an end of the parameter space (Self & Liang 1987). boundary : independence is on the boundary of the family. tau_k : the constrained (weighted) pseudo-maximum-likelihood τ̂. log_likelihood : Σ w log c at τ̂. n_obs, n_eff : positive-weight observations and Σ w.

submodel_lr_test

submodel_lr_test(full_family, sub_family, uv, weights=None) -> SubmodelLRTest

Likelihood-ratio test of a family's sub-model, one boundary parameter away.

H0: sub_family (the extra parameter at its sub-model value); H1: full_family. The (weighted) pseudo-likelihoods ℓ = Σ w log c(û, v̂; ·) are maximised separately over each model — the full family by the same joint MLE as :meth:CopulaVirt.fit(method='mle') on a multi-parameter family (:func:pmcprg.copulas._fit._fit_two_parameter_mle, which despite its name fits any number of extras) — and compared: LR = 2 max(ℓ_full − ℓ_sub, 0). The sub-model is fitted the same way, dispatched on how many parameters it has: the 1-D profile of :func:independence_lr_test when it is a plain one-parameter family (every nesting below except the last two), the same joint MLE as the full model when it is itself a two-parameter family — Tawn 1/2's own psi (Tawn-3 round, FR-9): profiling τ alone there would leave psi at the constructor's default instead of its own MLE, understating ℓ_sub and biasing the statistic upward.

Implemented nestings (module docstring, "Two-parameter sub-model tests", verifies each against the family's own module docstring): BB1 → Clayton, BB1 → Gumbel (CopulaGH), Student → Gauss (CopulaGaussian), Tawn type 1 or 2 → Gumbel (CopulaGH), BB6 → Joe, BB6 → Gumbel (CopulaGH), BB7 → Clayton, BB7 → Joe, BB8 → Joe, Tawn 3 → Tawn type 1 or 2. Every one is a boundary of the full family's admissible extra-parameter range, so LR → ½χ²₀ + ½χ²₁ (Self & Liang 1987), unlike the plain χ²₁ an interior sub-model would give.

Parameters:

Name Type Description Default
full_family the multi-parameter ``CopulaVirt`` subclass (or an
      instance) — H1.
required
sub_family
      it — H0, one or two parameters. Must be one of the
      nestings ``full_family`` registers above.
required
uv
required
weights
None

Raises:

Type Description
ValueError for a ``(full_family, sub_family)`` pair this module does not
implement, for a ``sub_family`` with more than two parameters, or when
either pseudo-likelihood is not finite at any evaluated parameter value.

SubmodelLRTest dataclass

SubmodelLRTest(family: str, submodel: str, statistic: float, p_value: float, null_distribution: str, boundary: bool, boundary_note: str, full_params: dict, sub_params: dict, log_likelihood_full: float, log_likelihood_sub: float, n_obs: int, n_eff: float)

Likelihood-ratio test of a one-parameter sub-model of a two-parameter family.

H0: the extra parameter sits at the sub-model's value (module docstring lists the three implemented nestings and why each is a boundary case); H1: the full two-parameter family.

Fields

family : class name of the two-parameter family (H1). submodel : class name of the one-parameter sub-model (H0). statistic : 2 (ℓ_full − ℓ_sub) ≥ 0, ℓ = Σ w log c at each model's own (weighted) pseudo-MLE. p_value : from null_distribution. null_distribution : '0.5*chi2(0) + 0.5*chi2(1)' for every pair this module implements — the sub-model sits at a boundary of the full family's admissible extra-parameter range in all three cases (Self & Liang 1987). boundary : always True here (see null_distribution). boundary_note : one sentence identifying which bound the sub-model is (module docstring). full_params : the full family's fitted {'tau_k': ..., extra: ...}. sub_params : the sub-model's fitted {'tau_k': ...}. log_likelihood_full, log_likelihood_sub : Σ w log c at each fit. n_obs, n_eff : positive-weight observations and Σ w.

mle_tau_discrepancy_test

mle_tau_discrepancy_test(family, uv, weights=None) -> MleTauDiscrepancyTest

Hausman-style test that the pseudo-MLE and Kendall's τ̂ agree (FR-7 b).

H0: the family is correctly specified and the sample uncontaminated, so both estimators are consistent for the same τ and τ̂_MLE − τ̂_τ = O_p(n^{-1/2}) is centred at 0. Kendall's τ̂ has a bounded influence function; the pseudo-likelihood score does not (Croux & Dehon 2010), so a few contaminating pairs move τ̂_MLE and leave τ̂_τ alone, and the standardised difference grows.

Both estimates are recomputed here from uv — the pseudo-MLE by the same bounded Brent search over the padded τ-range as :meth:CopulaVirt.fit(method='mle'), Kendall's τ̂ through the (weighted) empirical copula in the U-statistic form of :func:_tau_influence, which with unit weights and distinct pseudo-observations equals :func:scipy.stats.kendalltau exactly. The standard error of their difference comes from the two influence functions evaluated on the same sample, not from a difference of variances: see the module docstring, "MLE-vs-τ discrepancy diagnostic", for why the Hausman shortcut is not available here and what is computed instead.

Parameters:

Name Type Description Default
family
required
uv
required
weights optional (n,) non-negative frequency weights (module docstring,
  "Weights").
None

Returns:

Type Description
MleTauDiscrepancyTest — ``statistic`` is N(0, 1) under H0, ``p_value`` is
two-sided. At a boundary of the parameter space the statistic, its
standard error, the correlation and the p-value are NaN (``at_boundary``
set, ``boundary`` saying why): the pseudo-MLE is not asymptotically
normal there, so the contrast has no reference distribution.

Raises:

Type Description
ValueError for a family with no free parameter, an unusable

uv/weights, or a pseudo-likelihood that is not finite anywhere the search looked.

NotImplementedError for a two-parameter family: Kendall's τ does not

identify its second parameter, exactly as for :func:standard_errors with method='tau'.

MleTauDiscrepancyTest dataclass

MleTauDiscrepancyTest(family: str, statistic: float, p_value: float, difference: float, se_difference: float, tau_mle: float, tau_tau: float, theta_mle: float, theta_tau: float, se_tau_mle: float, se_tau_tau: float, correlation: float, at_boundary: bool, boundary: tuple, n_obs: int, n_eff: float)

Hausman-style contrast of the pseudo-MLE and the τ-inversion estimate.

H0: both estimators are consistent for the same τ — what correct specification of the family and uncontaminated data imply. See the module docstring ("MLE-vs-τ discrepancy diagnostic") for the construction of se_difference, which is not a difference of the two variances.

Fields

family : class name of the one-parameter family. statistic : (τ̂_MLE − τ̂_τ) / SE(τ̂_MLE − τ̂_τ); N(0, 1) under H0. p_value : two-sided, 2 Φ(−|statistic|). difference : τ̂_MLE − τ̂_τ. se_difference : its standard error, from the joint influence functions. tau_mle, tau_tau : the two estimates of τ. theta_mle, theta_tau : the family's native parameter at each — reported for information; the test is on the τ scale, which is equivalent to first order (module docstring). se_tau_mle, se_tau_tau : the two marginal standard errors of τ̂, identical to what :func:standard_errors returns for method='mle' and method='tau'. correlation : estimated Corr(τ̂_MLE, τ̂_τ) from the same influence functions — typically 0.95-0.99, which is exactly why the cross term may not be dropped. at_boundary : the pseudo-MLE sits on a boundary of the parameter space. statistic, p_value, se_difference and correlation are then NaN: neither the sandwich nor the normal limit applies there (Self & Liang 1987), the same convention :func:standard_errors follows. boundary : one sentence per boundary condition met (empty if interior). n_obs, n_eff : positive-weight observations and Σ w.

Robust fitting (density power divergence)

Weighted density-power-divergence estimation (audit FR-7), an alternative to MLE that trades a small efficiency loss on clean data for a bounded influence function under contamination.

dpd_fit

dpd_fit(family, uv, alpha: float, weights=None) -> DPDFit

Minimum-DPD estimate of a one-parameter family's τ (module docstring).

The objective :func:dpd_objective is minimised by the same bounded Brent search over the padded τ-range that :meth:CopulaVirt.fit(method='mle') uses, so alpha = 0 reproduces the pseudo-MLE to the optimiser's tolerance.

Parameters:

Name Type Description Default
family
required
uv
required
alpha float
  integral ``∫∫c^{1+α}`` becomes hard for tail-dependent families
  (module docstring); above 1 the call is refused.
required
weights optional (n,) non-negative frequency weights.
None

Returns:

Type Description
DPDFit

DPDFit dataclass

DPDFit(family: str, alpha: float, copula: object, tau_k: float, theta: float, objective: float, integral: float, integral_rel_error: float, variance: float, variance_rel_error: float, se_hint: float, converged: bool, n_obs: int, n_eff: float)

Minimum-density-power-divergence estimate of a one-parameter family.

Fields

family : class name. alpha : the tuning constant used (0 = pseudo-MLE). copula : the fitted instance (cls(tau_k=tau_k)). tau_k : τ̂_α. theta : the family's native parameter at τ̂_α. objective : H_n(τ̂_α; α) (module docstring) — the minimised value; comparable across τ at fixed α, not across α. integral : I(τ̂_α, α) = ∫∫ c^{1+α} at the estimate; 1.0 at α = 0. integral_rel_error : |I_refined − I_default| / I_refined at the estimate — the quadrature's own error estimate, 0.0 at α = 0 (where the integral is exactly 1) and 0.0 for the Gaussian copula (closed form). Compare it with se_hint before trusting the last digits of tau_k. variance : V̂_α = K_α/J_α², the known-margin asymptotic variance of √n_eff (τ̂_α − τ) under the model (module docstring, "Choosing α"); NaN when the two quadrature designs disagree, which is what a divergent ∫∫c^{1+2α} looks like (:func:_dpd_variance) — the estimate tau_k itself is unaffected. variance_rel_error : the relative disagreement behind that decision. se_hint : sqrt(variance / n_eff) — a hint, not a standard error: the rank-margin correction is not included (module docstring). converged : the bounded Brent search reported success and the objective is finite at the estimate. n_obs, n_eff : positive-weight observations and Σ w.

select_alpha

select_alpha(family, uv, weights=None, *, alphas=DPD_ALPHAS, pilot_alpha: float = 0.5) -> DPDAlphaSelection

Data-driven choice of α by the empirical-MSE criterion (module docstring).

Ĥ(α) = (τ̂_α − τ̂_P)² + V̂_α/n_eff is minimised over alphas, with the pilot τ̂_P the estimate at pilot_alpha. Read the module docstring's two caveats before quoting α̂: the criterion is Warwick & Jones (2005) with an explicit pilot rather than a reconstruction of Ghosh & Basu (2015)'s pilot-free variant, and V̂_α ignores the rank-margin correction.

Parameters:

Name Type Description Default
family
required
uv
required
weights
None
alphas
DPD_ALPHAS
pilot_alpha α of the pilot estimate (default 0.5).
0.5

Returns:

Type Description
DPDAlphaSelection

DPDAlphaSelection dataclass

DPDAlphaSelection(family: str, alpha: float, fit: DPDFit, alphas: tuple, tau_k: tuple, variance: tuple, criterion: tuple, pilot_alpha: float, pilot_tau: float, fits: dict, n_obs: int, n_eff: float)

Data-driven choice of the DPD tuning constant (module docstring).

Fields

family : class name. alpha : the selected α — the minimiser of criterion. fit : the :class:DPDFit at alpha. alphas : the grid searched, in the order of the arrays below. tau_k : τ̂_α on the grid. variance : V̂_α on the grid. criterion : (τ̂_α − τ̂_P)² + V̂_α/n_eff on the grid. pilot_alpha : the α of the pilot estimate τ̂_P. pilot_tau : τ̂_P itself. fits : every :class:DPDFit computed, keyed by α. n_obs, n_eff: positive-weight observations and Σ w.

DPD_ALPHAS module-attribute

DPD_ALPHAS: tuple[float, ...] = (0.0, 0.05, 0.1, 0.25, 0.5, 0.75)

Nonparametric comparison tool

Not a CopulaEnum family (no fit(), not selectable by ICE) — a nonparametric baseline to compare a fitted parametric family against.

EmpiricalBetaCopula

EmpiricalBetaCopula(data)

The empirical beta copula of Segers, Sibuya & Tsukahara (2017).

A standalone nonparametric comparison tool (AUDIT_COPULES FR-9) — not a :class:~pmcprg.copulas._base.CopulaVirt family: it has no parameter, is not in :class:~pmcprg.copulas._base.CopulaEnum, and is not used by ICE's family selection. See the module docstring for the full derivation, the tie convention, and the numerical-stability discussion.

Parameters:

Name Type Description Default
data (n, 2) array_like

Either raw bivariate observations or pseudo-observations already in [0, 1] — ranking is invariant to any strictly increasing per-column rescaling, so both give identical results. n ≥ 2 required (n ≥ 1 gives ranks a1 = a2 = 1, b1 = b2 = 1, i.e. the independence copula, which is degenerate as an estimate but not undefined; n ≥ 2 is required here so the estimator's own uniform-margin discrepancy — see the module docstring — is at least meaningful to quote).

required

Attributes:

Name Type Description
n int

Sample size.

a1, b1, a2, b2 (n,) ndarray

Beta shape parameters (R_i, n + 1 - R_i) for each coordinate; a1 + b1 == a2 + b2 == n + 1 exactly.

cdf_array

cdf_array(uv: ndarray) -> np.ndarray

C_n^β(u, v) at M point pairs — shape (M, 2) in, (M,) out.

cdf

cdf(u: float, v: float) -> float

Scalar convenience wrapper around :meth:cdf_array.

logpdf_array

logpdf_array(uv: ndarray) -> np.ndarray

log c_n^β(u, v) at M point pairs, via logsumexp over the n per-point terms (not floored; -inf where the true value underflows float64's range — see module docstring).

pdf_array

pdf_array(uv: ndarray) -> np.ndarray

c_n^β(u, v) = exp(logpdf_array(uv)) — see module docstring for why this is derived from the log form rather than the other way round (house convention, :class:CopulaVirt).

h1_array

h1_array(uv: ndarray) -> np.ndarray

h_1(v | u) = ∂C_n^β/∂u — the conditional cdf of V given U=u.

h2_array

h2_array(uv: ndarray) -> np.ndarray

h_2(u | v) = ∂C_n^β/∂v — the conditional cdf of U given V=v.

sample

sample(n: int = 500, seed: int | None = None) -> np.ndarray

Draw n samples from the mixture C_n^β exactly represents.

For each draw: pick a data index I ~ Uniform{0, …, n_data − 1}, then draw U ~ Beta(a1[I], b1[I]) and V ~ Beta(a2[I], b2[I]) independently (no Rosenblatt inversion needed — see module docstring).