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
Base classes¶
CopulaVirt ¶
CopulaVirt(class_name: str, params: dict)
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_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'
|
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:
The log-likelihood is the total Σ wᵢ log cᵢ and the BIC charges
|
None
|
pseudo_obs
|
bool
|
|
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
|
|
None
|
method
|
str
|
|
'tau'
|
weights
|
|
None
|
|
criterion
|
str
|
|
'aic'
|
pseudo_obs
|
passed to each ``cls.fit`` (see :meth:`fit`).
|
|
False
|
Returns:
| Type | Description |
|---|---|
FitBestResults
|
A |
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 overridesfit.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=Falsefor 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.
n_eff
property
¶
n_eff: float
Effective sample size: n_obs unweighted, Σw weighted (frequency weights).
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.
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.
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) whennis 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
+ 1is on the scale of the weights: they are frequency weights (an integer weight counts its observation that many times, the package's convention, as forn_eff = Σwand 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 withpseudo_obs=Truefor 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)
|
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
|
|
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
|
|
required |
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.
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
|
|
required | |
weights
|
optional (n,) non-negative frequency weights (module docstring).
|
|
None
|
method
|
str
|
|
'mle'
|
ranks
|
bool
|
|
True
|
Returns:
| Type | Description |
|---|---|
StandardErrors
|
|
Raises:
| Type | Description |
|---|---|
ValueError for a family without a free parameter, an unknown
|
method, |
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 = 1for 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
|
|
required |
sub_family
|
|
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,
|
|
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
|
|
NotImplementedError for a two-parameter family: Kendall's τ does not
|
identify its second parameter, exactly as for
:func: |
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
|
|
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.
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 |
cdf_array ¶
cdf_array(uv: ndarray) -> np.ndarray
C_n^β(u, v) at M point pairs — shape (M, 2) in, (M,) out.
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).