Extending pyvinecopulib: custom pair copulas and vines

The evaluators pyvinecopulib.core.Bicop / Vinecop and their PyTorch counterparts pyvinecopulib.torch.TorchBicop / TorchVinecop are concrete implementations of two array-agnostic contracts, BicopLike and VinecopLike. You can plug your own pair copula into a vine by implementing the contract — most easily by subclassing the canonical BicopBase / VinecopBase, which fill in almost everything from a few primitives.

For almost every use case a vine is simplified and unconditional: each pair copula has fixed parameters and sees only its two h-transformed arguments. That is the default, and where this notebook starts. It then shows how to select a structure from data with your own pairs, and the two advanced extensions VinecopBase also supports — a non-simplified vine (each pair also conditions on its edge’s conditioning set) and external covariates.

Everything here is plain NumPy; the same subclassing pattern works on PyTorch.

import numpy as np
from scipy.special import ndtr, ndtri  # standard-normal CDF / quantile

import pyvinecopulib as pv
from pyvinecopulib.core import BicopBase, NonSimplifiedContext, VinecopBase

rng = np.random.default_rng(0)

1. A custom pair copula

Subclass BicopBase and implement just three primitives — pdf, hfunc1, hfunc2 (the density and the two conditional CDFs). BicopBase supplies the rest: the inverse h-functions hinv1 / hinv2 (numerical), a simulate sampler, loglik, a density plot, and a __repr__.

The optional second argument x carries conditioning variables. We start with the ordinary case: with x=None the correlation is a fixed base_rho, so this is a plain Gaussian pair copula. (Section 3 switches x on to make it conditional.)

class GaussianBicop(BicopBase[np.ndarray]):
  """Gaussian pair copula; correlation optionally depends on covariates ``x``."""

  def __init__(self, *, base_rho=0.5, scale=0.6, rho_max=0.8):
    self._base_rho, self._scale, self._rho_max = base_rho, scale, rho_max

  def _rho(self, u, x):
    if x is None:  # ordinary (simplified) pair -> a fixed correlation
      return np.full(u.shape[0], self._base_rho)
    w = np.arange(1, x.shape[1] + 1)  # position weights (see section 3)
    return self._rho_max * np.tanh(self._scale * (x * w).sum(-1) / x.shape[1])

  def pdf(self, u, x=None):
    uc = np.clip(u, 1e-10, 1 - 1e-10)
    z1, z2 = ndtri(uc[:, 0]), ndtri(uc[:, 1])
    r = self._rho(u, x)
    om = 1 - r * r
    return np.exp(
      (2 * r * z1 * z2 - r * r * (z1**2 + z2**2)) / (2 * om)
    ) / np.sqrt(om)

  def hfunc1(self, u, x=None):  # P(U2 <= u2 | U1 = u1)
    uc = np.clip(u, 1e-10, 1 - 1e-10)
    z1, z2 = ndtri(uc[:, 0]), ndtri(uc[:, 1])
    r = self._rho(u, x)
    return ndtr((z2 - r * z1) / np.sqrt(1 - r * r))

  def hfunc2(self, u, x=None):  # P(U1 <= u1 | U2 = u2)
    uc = np.clip(u, 1e-10, 1 - 1e-10)
    z1, z2 = ndtri(uc[:, 0]), ndtri(uc[:, 1])
    r = self._rho(u, x)
    return ndtr((z1 - r * z2) / np.sqrt(1 - r * r))

  def flip(self):  # enables hosting in structure selection (section 2.1)
    # Gaussian is exchangeable: swapping the arguments leaves it unchanged.
    return GaussianBicop(
      base_rho=self._base_rho, scale=self._scale, rho_max=self._rho_max
    )

  def _simulate_uniform(self, n, qrng, seeds):  # enables simulate()
    return np.random.default_rng(seeds[0] if seeds else 0).uniform(size=(n, 2))

BicopBase fills in the rest from those three primitives — the inverse h-functions, a sampler, the log-likelihood, a density plot, and a repr:

cop = GaussianBicop(base_rho=0.6)
u = rng.uniform(0.05, 0.95, size=(5, 2))
print("pdf:     ", np.round(cop.pdf(u), 3))
print("hinv1:   ", np.round(cop.hinv1(u), 3), "(numerical inverse of hfunc1)")
print("loglik:  ", round(float(cop.loglik(u)), 3))
print("simulate:", np.round(cop.simulate(3, seeds=[1]), 3).tolist())
print(repr(cop))

pdf:      [0.952 2.689 1.68  1.281 0.91 ]
hinv1:    [0.402 0.021 0.915 0.719 0.852] (numerical inverse of hfunc1)
loglik:   1.612
simulate: [[0.512, 0.909], [0.144, 0.748], [0.312, 0.327]]
GaussianBicop()

2. Hosting it in a vine

VinecopBase implements the whole tree-by-tree cascade (pdf, rosenblatt, inverse_rosenblatt, simulate, cdf, loglik, plot, …) on top of a single required hook, _get_pair_copula, which returns the pair at (tree, edge). By default it assembles nothing extra per edge (SimplifiedContext), i.e. an ordinary unconditional vine.

class ListVinecop(VinecopBase[np.ndarray]):
  """A vine over a plain nested list of BicopLike pairs."""

  def __init__(self, pairs, structure, context=None):
    self._pairs = pairs
    self._bind_vine(structure, context)  # SimplifiedContext by default

  def _get_pair_copula(self, tree, edge):
    return self._pairs[tree][edge]

  def _simulate_uniform(self, n, qrng, seeds):
    return np.random.default_rng(seeds[0] if seeds else 0).uniform(
      size=(n, self.d)
    )


d = 4
structure = pv.RVineStructure.from_order([1, 2, 3, 4])
pairs = [
  [GaussianBicop(base_rho=r) for r in row]
  for row in ([0.5, 0.4, 0.3], [0.25, 0.2], [0.15])
]
vine = ListVinecop(pairs, structure)  # simplified + unconditional
print(vine, "| dim:", vine.dim, "| trees:", vine.trunc_lvl)

U = rng.uniform(0.05, 0.95, size=(1000, d))
print("log-likelihood:", round(float(vine.loglik(U)), 2))
sim = vine.simulate(500, seeds=[0])
print("simulate ->", sim.shape)

ListVinecop(dim=4, trunc_lvl=3, order=[1, 2, 3, 4]) | dim: 4 | trees: 3
log-likelihood: -99.97
simulate -> (500, 4)

A simplified vine of Gaussian pairs is exactly what the built-in pv.Vinecop builds, so hosting pv.Bicop Gaussian pairs in our VinecopBase subclass reproduces pv.Vinecop.from_structure to machine precision:

g = pv.families.gaussian
gauss = [
  [pv.Bicop(family=g, parameters=np.array([[r]])) for r in row]
  for row in ([0.5, 0.4, 0.3], [0.25, 0.2], [0.15])
]
ours = ListVinecop(gauss, structure)
ref = pv.Vinecop.from_structure(structure=structure, pair_copulas=gauss)
print("max |pdf difference|:", float(np.abs(ours.pdf(U) - ref.pdf(U)).max()))

max |pdf difference|: 1.1102230246251565e-15

2.1 Selecting a structure from data

So far the structure was given. VinecopBase.select builds one from data, running the same array-agnostic Dissmann selection pv.Vinecop uses: weight each candidate edge by Kendall’s τ, keep a maximum-dependence spanning tree, fit its pairs, and propagate their h-functions to the next tree. It returns the chosen RVineStructure together with the fitted pairs, reused — and reoriented via each pair’s flip — onto their finalized slots, so nothing is re-fit. You supply fit_edge(tree, edge, u_e, x_e), which fits one pair per edge; here it estimates a Gaussian correlation. (A pair hosted in selection therefore also needs a flip; evaluation along a fixed structure never calls it.)

def fit_edge(tree, edge, u_e, x_e):
  z = ndtri(np.clip(u_e, 1e-10, 1 - 1e-10))  # normal scores
  rho = float(np.corrcoef(z, rowvar=False)[0, 1])
  return GaussianBicop(base_rho=rho)


data = pv.to_pseudo_obs(
  rng.standard_normal((2000, d)) @ rng.standard_normal((d, d))
)
sel_structure, sel_pairs = VinecopBase.select(
  data, fit_edge, tree_criterion="tau"
)
selected = ListVinecop(sel_pairs, sel_structure)
print(
  "selected order:",
  list(sel_structure.order),
  "| trees:",
  sel_structure.trunc_lvl,
)
print("log-likelihood:", round(float(selected.loglik(data)), 1))

selected order: [2, 3, 1, 4] | trees: 3
log-likelihood: 1000.8

3. Advanced: non-simplified and conditional vines

The two extensions below cover the ~1% of cases that need them; a plain vine never has to touch them.

Non-simplified. Pass a NonSimplifiedContext: each pair copula c_{a,b;D} then also receives its edge’s conditioning-set values u_D (assembled by the cascade) as its x, so its correlation varies with the conditioning set. Our GaussianBicop already reads x in its _rho link, so the same kind of pair becomes a genuinely non-simplified vine just by switching the context. Because the conditioning variables are finalised before they are needed, inverse_rosenblatt still inverts rosenblatt exactly:

cond_pairs = [
  [GaussianBicop(scale=0.6) for _ in range(d - 1 - t)] for t in range(d - 1)
]
cond_vine = ListVinecop(cond_pairs, structure, NonSimplifiedContext())

W = cond_vine.rosenblatt(U)  # dependent -> independent uniforms
U_back = cond_vine.inverse_rosenblatt(W)  # and back
print("non-simplified round-trip max error:", float(np.abs(U_back - U).max()))
print("log-likelihood:", round(float(cond_vine.loglik(U)), 2))

non-simplified round-trip max error: 6.661338147750939e-15
log-likelihood: -223.9

External covariates. You can also pass an external covariate matrix x (row-aligned with u) to any evaluator; every pair sees it appended to its conditioning matrix. The density then shifts with the covariates:

X = rng.standard_normal(size=(1000, 2))
mean_at_zero = cond_vine.pdf(U, x=np.zeros((1000, 2))).mean()
mean_at_x = cond_vine.pdf(U, x=X).mean()
print("mean pdf, x = 0 :", round(float(mean_at_zero), 3))
print("mean pdf, x ~ N :", round(float(mean_at_x), 3))

mean pdf, x = 0 : 1.01
mean pdf, x ~ N : 1.087

What’s next

• The same subclassing pattern works on PyTorch: implement the primitives with torch ops (and make the class an nn.Module) for autograd and GPU.

• VinecopBase.select (section 2.1) chooses a simplified structure and fits its pairs in one pass. To fit pairs along a fixed structure — including a non-simplified one, edge by edge — drive VinecopBase.fit(structure, u, fit_edge, context=..., x=...).

• See pyvinecopulib.core.BicopBase / VinecopBase, BicopLike / VinecopLike, and SimplifiedContext / NonSimplifiedContext for the full contracts.