Localized Anomaly Detection via Differentiable D-vine Copulas

# Localized Anomaly Detection via Differentiable D-vine Copulas
Source: [https://arxiv.org/html/2607.25020](https://arxiv.org/html/2607.25020)
11institutetext:University of Trieste, Trieste, Italy
22institutetext:Idrostudi srl, Trieste, Italy
22email:nicholasandrea\.pearson@phd\.units\.itFrancesca Zanello Davide RussoLuca Bortolussi Francesca Cairoli

###### Abstract

Vine copulas provide a flexible framework for modeling complex multivariate distributions through a hierarchical decomposition into bivariate pair\-copulas\. Fitting a D\-vine requires selecting a copula family and parameter configuration for each pair\-copula from a set of candidates encoding different dependence patterns\. As the number of variables and candidate families increases, the number of possible configurations grows combinatorially\. Existing fitting procedures address this challenge through sequential greedy decisions, committing to a single locally optimal family at each step and potentially discarding configurations that would yield a better global fit\. To overcome this limitation, we propose a novel estimation framework that combines gradient\-based maximum likelihood estimation, enabled by our fully differentiable implementation, with a beam\-search strategy that maintains multiple competing D\-vine configurations throughout the fitting process\. This allows a broader exploration of the configuration space while remaining computationally tractable\. Building on the fitted D\-vine, we introduce a localized anomaly detection framework that exploits the hierarchical decomposition to produce both global anomaly scores and edge\-level explanations\. Statistical guarantees are provided through Mondrian conformal prediction, while the pair\-copula structure enables the localization of anomalies to specific variable relationships\. We evaluate the proposed framework on both benchmark and real\-world datasets, demonstrating its effectiveness for interpretable anomaly detection with uncertainty quantification\.

## 1Introduction

Anomaly Detection can be viewed as the task of identifying observations that deviate from the pattern exhibited by ordinary data\. It plays a prominent role in a wide range of real\-world applications such as fraud detection, industrial monitoring and healthcare\[chandola2009anomaly\]\. In these settings, the consequences of incorrect identification can be significant, making it essential not only to predict if an observation is anomalous, but also to quantify the reliability of the corresponding predictions\. It is often valuable to localize each anomaly, identifying the variables responsible for it and providing actionable insight into its underlying cause\. To address these needs, we approach the task at hand as a two\-stage problem\. First, we estimate the joint distribution of the observations under ordinary conditions to obtain a model of normal behavior through copula modeling\. Second, we build on this distribution, treating anomaly detection as a binary classification task, and employ a classifier whose predictions offer guarantees on their reliability through Conformal Prediction\[angelopoulos2023conformal\]\.

Copulas provide a flexible framework for modeling the joint distribution of a set of variables by separating marginal behavior from the dependence structure\[sklar1959fonctions\]\. Standard multivariate copula scale poorly to high\-dimensions\[nelsen2006introduction\], a limitation overcome by pair\-copula constructions\[aas2009pair\], which decompose the multivariate distributions into a hierarchy of bivariate copulas\. Out of these, we focus on the D\-vine, which arranges the variables along a path ordering\[czado2019analyzing\]\. Fitting a D\-vine requires selecting a family and parameters for each pair\-copula, a combinatorial problem that existing methods address either by exploring too few configurations\[dissmann2013selecting\]or at high computational cost\[chapon2023imputation\]\.

We propose a methodology that combines an efficient D\-vine fitting procedure with a conformal anomaly detection scheme providing per\-class reliability guarantees\. Our contributions are:

1. i\.A fully differentiable PyTorch implementation of bivariate copula families, supporting gradient\-based parameter estimation and GPU acceleration\.
2. ii\.A beam search procedure for D\-vine fitting that, at each pair\-copula, evaluates all candidate families and propagates forward only the most promising configurations\. Several candidate paths are evaluated simultaneously, allowing for a broader exploration of the configuration space without incurring the cost of exhaustive search\.
3. iii\.An anomaly detection procedure framing detection as a classification task between ordinary and anomalous behavior\. Mondrian conformal prediction is applied to obtain prediction regions that are guaranteed to contain the true label with a prescribed confidence level for each class\. The hierarchical structure of D\-vines is then exploited to localize each anomaly to individual pair\-copulas in the decomposition and thus to the responsible variables\.

Related Works\.Several approaches have been proposed for*fitting D\-vines*, with the standard approach of\[dissmann2013selecting\]selecting a single configuration in a greedy procedure which limits the number of explored configurations\. Bayesian alternatives like the RJMCMC procedure of\[chapon2023imputation\]fully explore the configuration space but suffer from high computational costs and convergence issues\. Examples of*copula based anomaly detection*include CoCAI\[pearson2025cocai\]which combines copula modeling with conformal prediction for time series, and COPOD\[li2020copod\]which uses empirical copulas for outlier detection\. Lastly,\[horvath2020copula\]investigates anomalous dependence patterns through unconditional pair\-copulas defined on a spanning\-tree structure\.

## 2Background

### 2\.1Copula Theory

LetX=\(X1,…,Xd\)X=\(X\_\{1\},\\dots,X\_\{d\}\)be a random vector with joint cumulative distribution function \(CDF\)FFand continuous marginal CDFsF1,…,FdF\_\{1\},\\dots,F\_\{d\}\. By the probability integral transform, the variablesUi=Fi​\(Xi\)U\_\{i\}=F\_\{i\}\(X\_\{i\}\)withi=1,…,d,i=1,\\dots,d,are uniformly distributed on\[0,1\]\[0,1\]\[yan2007enjoy\]\. The joint distribution of\(U1,…,Ud\)\(U\_\{1\},\\dots,U\_\{d\}\)is called a*copula and characterizes the dependence structure ofXXindependently of its marginals*\[nelsen2006introduction\]\. The central result of copula theory is Sklar’s theorem\[sklar1959fonctions\]\.

###### Theorem 2\.1\(Sklar’s Theorem\)

LetFFbe add\-dimensional CDF with continuous marginalsF1,…,FdF\_\{1\},\\dots,F\_\{d\}\. Then there exists a unique copulaC:\[0,1\]d→\[0,1\]C:\[0,1\]^\{d\}\\to\[0,1\]such that

F​\(x1,…,xd\)=C​\(F1​\(x1\),…,Fd​\(xd\)\)\.F\(x\_\{1\},\\dots,x\_\{d\}\)=C\\left\(F\_\{1\}\(x\_\{1\}\),\\dots,F\_\{d\}\(x\_\{d\}\)\\right\)\.\(1\)Conversely, for any copulaCCand marginal CDFsF1,…,FdF\_\{1\},\\dots,F\_\{d\}, Equation \([1](https://arxiv.org/html/2607.25020#S2.E1)\) defines a validdd\-dimensional distribution\.

Sklar’s theorem separates marginal behavior from dependence\. IfCCis absolutely continuous, its copula density is defined as

c​\(u1,…,ud\)=∂dC​\(u1,…,ud\)∂u1​⋯​∂ud\.c\(u\_\{1\},\\dots,u\_\{d\}\)=\\frac\{\\partial^\{d\}C\(u\_\{1\},\\dots,u\_\{d\}\)\}\{\\partial u\_\{1\}\\cdots\\partial u\_\{d\}\}\.\(2\)If the marginals are absolutely continuous with densitiesf1,…,fdf\_\{1\},\\dots,f\_\{d\}, the joint density factorizes as

f​\(x1,…,xd\)=c​\(F1​\(x1\),…,Fd​\(xd\)\)​∏i=1dfi​\(xi\)\.f\(x\_\{1\},\\dots,x\_\{d\}\)=c\\left\(F\_\{1\}\(x\_\{1\}\),\\dots,F\_\{d\}\(x\_\{d\}\)\\right\)\\prod\_\{i=1\}^\{d\}f\_\{i\}\(x\_\{i\}\)\.\(3\)This decomposition forms the basis of likelihood\-based inference and, in particular, of pair\-copula constructions\[aas2009pair\]\.

#### 2\.1\.1Pair\-Copula Constructions\.

Direct modeling with high\-dimensional copulas is often impractical due to limited flexibility and restrictive dependence assumptions\[min2010bayesian\]\. Pair\-copula constructions \(PCCs\) address this by decomposing a multivariate copula into a collection of bivariate copulas, enabling flexible and scalable modeling of complex dependence structures\[aas2009pair\]\. Bivariate copulas model the dependence between two variables, and are also referred to as pair\-copulas\. They are typically classified into the elliptical and Archimedean families\. Elliptical copulas, such as the Gaussian and Student\-t copulas, are characterized by symmetric dependence patterns\. In contrast, Archimedean copulas, including the Gumbel and Clayton copulas, capture asymmetric and tail\-dependent behaviors\[dissmann2013selecting\]\. The dependence captured by a pair\-copula is characterized by a set of family\-specific parameters governing the strength of the relationship between the corresponding variables\. The probability that the two variables simultaneously take extremely high or low values can be summarized by the copula’s tail dependence coefficients, admitting values in\[0,1\]\[0,1\], with values close to zero indicating tail independence and those close to one indicating strong dependence\[nelsen2006introduction\]\. A more detailed characterization of tail dependence for each copula family is provided in AppendixLABEL:app:families\. A key property of bivariate copulas is that conditional distributions can be expressed through the partial derivatives of the copula\[schepsmeier2014derivatives\]\. For notational convenience, we denote these conditional distribution functions as h\-functions, with:

h1C​\(u1,u2\)=∂C​\(u1,u2\)∂u2,h2C​\(u1,u2\)=∂C​\(u1,u2\)∂u1\.h\_\{1\}^\{C\}\(\{u\_\{1\}\},\{u\_\{2\}\}\)=\\frac\{\\partial C\(u\_\{1\},u\_\{2\}\)\}\{\\partial u\_\{2\}\},\\qquad h\_\{2\}^\{C\}\(\{u\_\{1\}\},\{u\_\{2\}\}\)=\\frac\{\\partial C\(u\_\{1\},u\_\{2\}\)\}\{\\partial u\_\{1\}\}\.These functions play a central role in PCCs, enabling the recursive construction of higher\-dimensional dependence models starting from bivariate copulas\.

![Refer to caption](https://arxiv.org/html/2607.25020v1/Figures/DiagrammaDVine.png)Figure 1:D\-vine decomposition ofd=5d=5variables\. Each row is a treeT1T\_\{1\}\-T4T\_\{4\}and each edge a pair\-copula\. Nodes inT1T\_\{1\}are the modeled variables, while in subsequent trees they correspond to the edges of the previous tree\. Dashed arrows represent the h\-functions \(or pseudo observationsuiu\_\{i\}inT1T\_\{1\}\) that propagate data from one tree to the next\.A PCC indddimensions consists ofd​\(d−1\)/2d\(d\-1\)/2pair\-copulas\[aas2009pair\]\. A common PCC structure is the*D\-vine*, which arranges these pair\-copulas as edges acrossd−1d\-1linked trees \(Fig\.[1](https://arxiv.org/html/2607.25020#S2.F1)\)\. The first tree contains unconditional pair\-copulas between adjacent variables in a chosen ordering\. Higher\-level trees introduce conditional pair\-copulas between variables that are progressively farther apart, conditioned on the variables lying between them\. The required conditional distributions are computed recursively through h\-functions\[hobaek2013parameter\]\.

This construction yields the factorization of add\-dimensional density as:

f​\(x1,…,xd\)=∏i=1dfi​\(xi\)⋅∏ℓ=1d−1∏i=1d−ℓci,i\+ℓ∣Di,ℓ​\(Fi∣Di,ℓ,Fi\+ℓ∣Di,ℓ\)f\(x\_\{1\},\\dots,x\_\{d\}\)=\\prod\_\{i=1\}^\{d\}f\_\{i\}\(x\_\{i\}\)\\cdot\\prod\_\{\\ell=1\}^\{d\-1\}\\prod\_\{i=1\}^\{d\-\\ell\}c\_\{i,i\+\\ell\\mid D\_\{i,\\ell\}\}\(F\_\{i\\mid D\_\{i,\\ell\}\},F\_\{i\+\\ell\\mid D\_\{i,\\ell\}\}\)\(4\)whereDi,ℓ=i\+1,…,i\+ℓ−1D\_\{i,\\ell\}=\{i\+1,\\dots,i\+\\ell\-1\}denotes the set of intermediate variables conditioning the pair\-copulaci,i\+ℓ∣Di,ℓc\_\{i,i\+\\ell\\mid D\_\{i,\\ell\}\}\[czado2019analyzing\]\. For example, whend=3d=3:

f​\(x1,x2,x3\)=f1​\(x1\)​f2​\(x2\)​f3​\(x3\)​c12​\(u1,u2\)​c23​\(u2,u3\)​c13\|2​\(u1\|2,u3\|2\),f\(x\_\{1\},x\_\{2\},x\_\{3\}\)=f\_\{1\}\(x\_\{1\}\)f\_\{2\}\(x\_\{2\}\)f\_\{3\}\(x\_\{3\}\)\\,c\_\{12\}\(u\_\{1\},u\_\{2\}\)\\,c\_\{23\}\(u\_\{2\},u\_\{3\}\)\\,c\_\{13\|2\}\(u\_\{1\|2\},u\_\{3\|2\}\),whereu1\|2u\_\{1\|2\}andu3\|2u\_\{3\|2\}are obtained via h\-functions\.

Fitting a D\-vine requires selecting both a copula family and its parameters for each of thed​\(d−1\)/2d\(d\-1\)/2pair\-copulas\[dissmann2013selecting\]\. Since the number of possible configurations grows combinatorially with the number of candidate families and variable orderings, exhaustive search quickly becomes unfeasible\[chapon2023imputation\]\. Moreover, pair\-copulas cannot be optimized independently: the conditional distributions used in higher\-level trees depend on the h\-functions of lower\-level trees, causing early modeling decisions to propagate throughout the vine\[hobaek2013parameter\]\.

### 2\.2Conformal Predictions

Conformal prediction provides distribution\-free uncertainty quantification by associating each test sample𝐱\\mathbf\{x\}with a prediction region rather than a single point prediction\[angelopoulos2023conformal\]\. Given a trained classifier, a separate calibration set𝒟cal\\mathcal\{D\}\_\{\\mathrm\{cal\}\}is used to compute non\-conformity scores𝚗𝚌𝚏i\\mathtt\{ncf\}\_\{i\}for calibration examples\. Given a significance levelα\\alpha, for a candidate labelyy, the test point\(x,y\)\(x,y\)is assigned a non\-conformity score𝚗𝚌𝚏​\(x,y\)\\mathtt\{ncf\}\(x,y\), which is compared against the threshold computed on the calibration scores\[fontana2023conformal\]:

τα=⌈\(∣𝒟cal∣\+1\)​\(1−α\)⌉∣𝒟cal∣​\-th quantile of​\{𝚗𝚌𝚏i:i∈𝒟cal\}\.\\tau\_\{\\alpha\}=\\frac\{\\lceil\(\\mid\\mathcal\{D\}\_\{\\mathrm\{cal\}\}\\mid\+1\)\(1\-\\alpha\)\\rceil\}\{\\mid\\mathcal\{D\}\_\{\\mathrm\{cal\}\}\\mid\}\\text\{\-th quantile of \}\\\{\\mathtt\{ncf\}\_\{i\}:i\\in\\mathcal\{D\}\_\{\\mathrm\{cal\}\}\\\}\.The prediction region is then defined as:Γα​\(x\)=\{y∈\{0,1\}:𝚗𝚌𝚏​\(x,y\)≤τα\}\.\\Gamma\_\{\\alpha\}\(x\)=\\\{y\\in\\\{0,1\\\}:\\mathtt\{ncf\}\(x,y\)\\leq\\tau\_\{\\alpha\}\\\}\.

In binary classification,Γα​\(x\)\\Gamma\_\{\\alpha\}\(x\)may contain one label \(confident prediction\), both labels \(uncertainty\), or, more rarely, no label\.*Mondrian*conformal prediction computes calibration scores separately for each class, yielding class\-conditional coverage guarantees\[vovk2005algorithmic\]\. Under the assumption of exchangeability between calibration and test observations belonging to the same class:

ℙ\(x,y\)​\(y∈Γα​\(x\)\)≥1−α,y∈\{0,1\},\\mathbb\{P\}\_\{\(x,y\)\}\\bigl\(y\\in\\Gamma\_\{\\alpha\}\(x\)\\bigr\)\\geq 1\-\\alpha,\\qquad y\\in\\\{0,1\\\},meaning that the prediction set contains the true label with probability at least1−α1\-\\alphafor each class individually\. These guarantees hold in finite samples and do not rely on assumptions about the underlying data distribution\[fontana2023conformal\]\.

## 3Methodology

This section presents our differentiable beam\-search procedure for fitting D\-vine copulas\. We address two key challenges of D\-vine model selection: the combinatorial growth of possible family assignments and the error propagation induced by greedy sequential fitting\. The method combines gradient\-based maximum likelihood estimation of pair\-copula parameters with a beam\-search strategy that maintains multiple candidate D\-vine configurations during the fitting process\.

### 3\.1Fitting D\-Vine Copulas

We assume access to a dataset of continuous variables suitable for copula modeling, where each observation is labeled as either ordinary or anomalous\. The data are partitioned into three disjoint subsets: training, calibration and test\. We denote with𝐗ry\\mathbf\{X\}\_\{r\}^\{y\}the subset corresponding to splitr∈\{tr,ca,te\}r\\in\\\{\\mathrm\{tr\},\\mathrm\{ca\},\\mathrm\{te\}\\\}and labely⊂\{0,1\}y\\subset\\\{0,1\\\}, where0denotes ordinary observations and11anomalous ones\. Thus,𝐗tr\(0,1\)\\mathbf\{X\}\_\{\\mathrm\{tr\}\}^\{\(0,1\)\}are used for model fitting,𝐗ca\(0,1\)\\mathbf\{X\}\_\{\\mathrm\{ca\}\}^\{\(0,1\)\}for calibration, and𝐗te\(0,1\)\\mathbf\{X\}\_\{\\mathrm\{te\}\}^\{\(0,1\)\}for evaluation\.

Given a dataset𝐗≡𝐗tr0∈ℝn×d\\mathbf\{X\}\\equiv\\mathbf\{X\}^\{0\}\_\{\\mathrm\{tr\}\}\\in\\mathbb\{R\}^\{n\\times d\}ofnn*ordinary*observations of dimensiondd, we first transform each variable to the unit interval\. For each variable𝐗⋅j\\mathbf\{X\}\_\{\\cdot j\}, we estimate its marginal CDFF^j\\hat\{F\}\_\{j\}using Gaussian kernel density estimation \(KDE\)\[silverman2018density\]and obtain pseudo\-observations𝐔i​j=F^j​\(𝐗i​j\)\\mathbf\{U\}\_\{ij\}=\\hat\{F\}\_\{j\}\(\\mathbf\{X\}\_\{ij\}\)\. By the probability integral transform\[yan2007enjoy\], the resulting matrix𝐔∈\(0,1\)n×d\\mathbf\{U\}\\in\(0,1\)^\{n\\times d\}contains approximately uniform marginals and forms the input to the D\-vine fitting procedure\.

The objective is to fit a D\-vine copula by selecting an appropriate pair\-copula family for each edge and estimating the corresponding parameters\. At a high level, the procedure alternates between local fitting and global search\. For every edge in the current tree, all candidate copula families are fitted independently by gradient\-based maximum likelihood estimation\. Rather than committing immediately to the best\-fitting family, statistically indistinguishable alternatives are retained\. These alternatives generate multiple extensions of the current partial D\-vine configurations\. A beam\-search strategy then preserves only the most promising configurations according to their cumulative log\-likelihood\. Once all trees have been processed, the surviving complete D\-vine configurations undergo a final joint parameter refinement\.

We consider a candidate family setℱ=\{F1,…,FK\}\\mathcal\{F\}=\\\{F\_\{1\},\\dots,F\_\{K\}\\\}\. For each pair\-copula in the D\-vine, the goal is to select a familyFk∈ℱF\_\{k\}\\in\\mathcal\{F\}and estimate its associated parameter vectorθk\\theta\_\{k\}\. The inputs to a generic pair\-copula are denoted byu,v∈\[0,1\]nu,v\\in\[0,1\]^\{n\}\. In the first tree, these correspond to pseudo\-observations obtained directly from the data, whereas in higher\-order trees, they are produced through the cascade of h\-functions\. A generic pair\-copula is denoted byCi​j∣DC\_\{ij\\mid D\}, whereiiandjjidentify the coupled variables andDDis the conditioning set\. To simplify notation, each pair\-copula is uniquely identified by its edge position\(t,p\)\(t,p\)in the D\-vine, wheret∈\{1,…,d−1\}t\\in\\\{1,\\dots,d\-1\\\}is the tree index andp∈\{1,…,d−t\}p\\in\\\{1,\\dots,d\-t\\\}is the edge position within treeTtT\_\{t\}\. We use the shorthand notationc\(t,p\)c\_\{\(t,p\)\},h1\(t,p\)h\_\{1\}^\{\(t,p\)\}andh2\(t,p\)h\_\{2\}^\{\(t,p\)\}for the corresponding copula density and h\-functions associated with the underlying pair\-copulaCi​j∣DC\_\{ij\\mid D\}\(see Fig\.[2](https://arxiv.org/html/2607.25020#S3.F2)\)\.

Since the D\-vine is fitted sequentially across treesT1,…,Td−1T\_\{1\},\\dots,T\_\{d\-1\}, we define the set of edges processed up to treeTtT\_\{t\}:

ℰ→t:=\{\(t′,p′\)∣t′=1,…,t,;p′=1,…,d−t′\}\.\\mathcal\{E\}\_\{\\to t\}:=\\\{\(t^\{\\prime\},p^\{\\prime\}\)\\mid t^\{\\prime\}=1,\\dots,t,;p^\{\\prime\}=1,\\dots,d\-t^\{\\prime\}\\\}\.
A*beam state*BtiB\_\{t\}^\{i\}represents a partial D\-vine configuration after processing treeTtT\_\{t\}\. It stores the selected family and parameter values for all edges inℰ→t\\mathcal\{E\}\_\{\\to t\}:

Bti:=\{b\(t′,p′\)∣\(t′,p′\)∈ℰ→t\},B\_\{t\}^\{i\}:=\\big\\\{b\_\{\(t^\{\\prime\},p^\{\\prime\}\)\}\\mid\(t^\{\\prime\},p^\{\\prime\}\)\\in\\mathcal\{E\}\_\{\\to t\}\\big\\\},whereb\(t′,p′\)=\(Fk,θk\)b\_\{\(t^\{\\prime\},p^\{\\prime\}\)\}=\(F\_\{k\},\\theta\_\{k\}\)is the selected family–parameter configuration for edge\(t′,p′\)\(t^\{\\prime\},p^\{\\prime\}\)\. Note that every edge\(t′,p′\)\(t^\{\\prime\},p^\{\\prime\}\)may be assigned to a different familyFkF\_\{k\}\. The setℬt=\{Bt1,Bt2,…\}\\mathcal\{B\}\_\{t\}=\\\{B\_\{t\}^\{1\},B\_\{t\}^\{2\},\\dots\\\}contains all active beam states after processing treeTtT\_\{t\}\. For edge\(t,p\)\(t,p\)and family–parameter configurationb\(t,p\)=\(Fk,θk\)b\_\{\(t,p\)\}=\(F\_\{k\},\\theta\_\{k\}\), for somek∈\{1,…,K\}k\\in\\\{1,\\dots,K\\\}, the*pointwise log\-likelihood*contribution of observationmmis:

ℓ​\(b\(t,p\);um,vm\)=log⁡\(c\(t,p\)​\(um,vm;θk\)\),\\ell\(b\_\{\(t,p\)\};u\_\{m\},v\_\{m\}\)=\\log\\big\(c\_\{\(t,p\)\}\(u\_\{m\},v\_\{m\};\\theta\_\{k\}\)\\big\),\(5\)wherec\(t,p\)​\(⋅,⋅;θk\)c\_\{\(t,p\)\}\(\\cdot,\\cdot;\\theta\_\{k\}\)is the copula density of familyFkF\_\{k\}\. The*edge\-specific log\-likelihood*is then computed as:

ℓ​\(b\(t,p\);𝐔\)=∑m=1nℓ​\(b\(t,p\);um,vm\)\.\\ell\(b\_\{\(t,p\)\};\\mathbf\{U\}\)=\\sum\_\{m=1\}^\{n\}\\ell\(b\_\{\(t,p\)\};u\_\{m\},v\_\{m\}\)\.\(6\)The*score*of a beam state is defined as the cumulative log\-likelihood of all edges it contains:

ℓ​\(Bti;𝐔\)=∑\(t′,p′\)∈ℰ→tℓ​\(b\(t′,p′\);𝐔\)\.\\ell\(B\_\{t\}^\{i\};\\mathbf\{U\}\)=\\sum\_\{\(t^\{\\prime\},p^\{\\prime\}\)\\in\\mathcal\{E\}\_\{\\to t\}\}\\ell\\left\(b\_\{\(t^\{\\prime\},p^\{\\prime\}\)\};\\mathbf\{U\}\\right\)\.\(7\)
The search breadth is controlled by two user\-defined hyperparameters: the maximum number of beam states retained after each treeβmax\\beta\_\{\\max\}, and the maximum number of candidate configurations retained for each edgeωmax\\omega\_\{\\max\}\.

The algorithm is initialized with a single empty beam stateB0∅B\_\{0\}^\{\\emptyset\}, withℬ0=\{B0∅\}\\mathcal\{B\}\_\{0\}=\\\{B\_\{0\}^\{\\emptyset\}\\\}\. Fitting proceeds sequentially over the treesT1,…,Td−1T\_\{1\},\\dots,T\_\{d\-1\}\. At iterationtt, each active beam stateBt−1i∈ℬt−1B\_\{t\-1\}^\{i\}\\in\\mathcal\{B\}\_\{t\-1\}is expanded independently\. For every edge\(t,p\)\(t,p\)in the current tree and for every candidate familyFk∈ℱF\_\{k\}\\in\\mathcal\{F\}, parameters are estimated through gradient\-based maximization of the edge\-specific log\-likelihood defined in \([6](https://arxiv.org/html/2607.25020#S3.E6)\):

θ^k=arg⁡maxθk⁡ℓ​\(b\(t,p\);𝐔\)\.\\hat\{\\theta\}\_\{k\}=\\arg\\max\_\{\\theta\_\{k\}\}\\ \\ell\(b\_\{\(t,p\)\};\\mathbf\{U\}\)\.\(8\)The fitted configurationb^\(t,p\)=\(Fk,θ^k\)\\hat\{b\}\_\{\(t,p\)\}=\(F\_\{k\},\\hat\{\\theta\}\_\{k\}\)is associated with edge log\-likelihoodℓ​\(b^\(t,p\);𝐔\)\.\\ell\(\\hat\{b\}\_\{\(t,p\)\};\\mathbf\{U\}\)\.After fitting all theKKcandidate families, the configuration with the highest edge log\-likelihood is identified as the leading candidate\. However, selecting only the top\-scoring family may discard alternatives that are statistically indistinguishable and could lead to better overall D\-vine configurations\.

To detect such cases, we compare the leading candidate against the other configurations using either the Vuong test\[vuong1989likelihood\]or the Clarke test\[clarke2007simple\]\. Both tests assess whether differences in likelihood are statistically significant\. If the leading candidate is significantly better than every alternative, the edge commits to it\. Otherwise, all non\-rejected alternatives are retained\. Since the Clarke test is generally more conservative, it typically produces larger candidate sets\.

To limit combinatorial growth, at mostωmax\\omega\_\{\\max\}candidate configurations are retained for each edge, prioritizing those with the highest edge log\-likelihood\. This produces the candidate set

𝒜\(t,p\)i=\{b^\(t,p\)1,…,b^\(t,p\)ω\},ω≤ωmax,\\mathcal\{A\}\_\{\(t,p\)\}^\{i\}=\\\{\\hat\{b\}\_\{\(t,p\)\}^\{1\},\\dots,\\hat\{b\}\_\{\(t,p\)\}^\{\\omega\}\\\},\\qquad\\omega\\leq\\omega\_\{\\max\},for edge\(t,p\)\(t,p\)under parent beam stateBt−1iB\_\{t\-1\}^\{i\}\. Once candidate sets have been obtained for all edges in the current tree, they are combined through the Cartesian product𝒜ti=𝒜\(t,1\)i×⋯×𝒜\(t,d−t\)i\.\\mathcal\{A\}\_\{t\}^\{i\}=\\mathcal\{A\}\_\{\(t,1\)\}^\{i\}\\times\\dots\\times\\mathcal\{A\}\_\{\(t,d\-t\)\}^\{i\}\.Each elementa∈𝒜tia\\in\\mathcal\{A\}\_\{t\}^\{i\}defines a possible extension of the parent beam stateBt−1iB\_\{t\-1\}^\{i\}\. The resulting beam score is obtained by adding the cumulative score of the parent beam to the edge log\-likelihoods contributed by the configurations inaa\. Repeating this process for every active beam state inℬt−1\\mathcal\{B\}\_\{t\-1\}generates a pool of candidate beam states for treeTtT\_\{t\}\. From this pool, only theβmax\\beta\_\{\\max\}states with the highest beam score are retained, forming the beam setℬt\\mathcal\{B\}\_\{t\}passed to the next iteration\.

After processing all trees, the final beam setℬd−1\\mathcal\{B\}\_\{d\-1\}contains complete D\-vine configurations\. Although these configurations achieve high cumulative log\-likelihood, pair\-copula parameters have been estimated locally while treating the outputs of previously fitted copulas as fixed\. Since the inputs of a pair\-copula in treeTtT\_\{t\}are generated through the h\-function cascade of preceding trees, parameter updates can affect all downstream pair\-copulas\.

To account for these dependencies, each beam state inℬd−1\\mathcal\{B\}\_\{d\-1\}undergoes a final*joint refinement*phase\. Keeping the selected copula families fixed, all pair\-copula parameters are optimized simultaneously by maximizing the full D\-vine log\-likelihood\. Optimization is initialized from the locally fitted parameters and performed end\-to\-end, allowing gradients to propagate through the complete h\-function cascade and explicitly accounting for cross\-tree parameter interactions\.

The beam state with the highest refined log\-likelihood is selected as the final fitted D\-vine configuration:

Bd−1∗=arg⁡maxB∈ℬd−1⁡ℓ​\(B;𝐔\)\.B\_\{d\-1\}^\{\*\}=\\arg\\max\_\{B\\in\\mathcal\{B\}\_\{d\-1\}\}\\ell\(B;\\mathbf\{U\}\)\.\(9\)Fig\.[2](https://arxiv.org/html/2607.25020#S3.F2)illustrates the full fitting procedure on an example withd=3d=3variables\.

##### Variable Ordering\.

The ordering of variables determines the unconditional pair\-copulas in the first treeT1T\_\{1\}and therefore strongly influences the quality of the resulting D\-vine\. Following the approach of\[dissmann2013selecting\], we construct the ordering by placing strongly dependent variables adjacent in the sequence, quantifying dependence between each pair of variables through the empirical Kendall’s Tau\[nelsen2006introduction\]\. The final ordering is the path through theddvariables maximizing the sum of absolute pairwise Kendall’s Tau over adjacent pairs\.

![Refer to caption](https://arxiv.org/html/2607.25020v1/Figures/DiagrammaBeam.png)Figure 2:Graphical representation of the proposed D\-vine fitting methodology ford=3d=3variables, withβmax=2\\beta\_\{\\text\{max\}\}=2andωmax=2\\omega\_\{\\text\{max\}\}=2\. Trees are processed sequentially from an empty beam stateB0∅B\_\{0\}^\{\\emptyset\}\. Each edge inT1T\_\{1\}is individually fitted, keeping at most the bestωmax\\omega\_\{\\text\{max\}\}configurations for each node\. These are combined via their Cartesian product into candidate beam statesB11,…,B14B\_\{1\}^\{1\},\\dots,B\_\{1\}^\{4\}, of which the bestβmax\\beta\_\{\\text\{max\}\}formℬ1\\mathcal\{B\}\_\{1\}\. The procedure repeats atT2T\_\{2\}, with the best final beams undergoing joint refinement yielding the final configurationB2∗B\_\{2\}^\{\*\}\. Crossed out states are discarded configurations\.
##### Anomaly\-aware fitting\.

The D\-vine is fitted on𝐗tr0\\mathbf\{X\}\_\{\\mathrm\{tr\}\}^\{0\}using the procedure described above\. When anomalous training observations𝐗tr1\\mathbf\{X\}\_\{\\mathrm\{tr\}\}^\{1\}are available, we augment the fitting objective with a penalization term that discourages high likelihood values on anomalous samples\. Specifically, Eq\.[6](https://arxiv.org/html/2607.25020#S3.E6)becomes

ℓpen​\(b\(t,p\);𝐔\)=ℓ​\(b\(t,p\);𝐔0\)−λ​ℓ​\(b\(t,p\);𝐔1\),\\ell^\{\\mathrm\{pen\}\}\(b\_\{\(t,p\)\};\\mathbf\{U\}\)=\\ell\(b\_\{\(t,p\)\};\\mathbf\{U\}^\{0\}\)\-\\lambda\\,\\ell\(b\_\{\(t,p\)\};\\mathbf\{U\}^\{1\}\),\(10\)where𝐔0\\mathbf\{U\}^\{0\}and𝐔1\\mathbf\{U\}^\{1\}are the pseudo\-observations derived from𝐗tr0\\mathbf\{X\}\_\{\\mathrm\{tr\}\}^\{0\}and𝐗tr1\\mathbf\{X\}\_\{\\mathrm\{tr\}\}^\{1\}, respectively, andλ≥0\\lambda\\geq 0controls the penalization strength\. This encourages the fitted D\-vine to capture dependence structures characteristic of ordinary observations while reducing its fit to anomalous ones\.

### 3\.2Anomaly Detection

We now describe how the fitted D\-vine is used to perform anomaly detection\. The proposed procedure leverages the hierarchical structure of the vine decomposition to produce both a global anomaly score and per\-edge scores that facilitate the identification of the variable relationships responsible for a detected anomaly\.

##### Edge\-level anomaly scores\.

The anomaly detection criterion is based on the negative log\-copula density associated with each pair\-copula in the fitted D\-vine under configurationBd−1∗B\_\{d\-1\}^\{\*\}\. Given an observation𝐮m\\mathbf\{u\}\_\{m\}, we compute the edge\-wise scores−ℓ​\(b\(t,p\)∗;𝐮m\)\-\\ell\(b^\{\*\}\_\{\(t,p\)\};\\mathbf\{u\}\_\{m\}\)\. These values are large when the observed dependence pattern deviates from that expected under the fitted pair\-copula and small when the two are consistent\. To establish a baseline under ordinary conditions, we evaluate these scores on the ordinary training set𝐗𝚝𝚛0\\mathbf\{X\}\_\{\\mathtt\{tr\}\}^\{0\}\. Since pair\-copulas belonging to different families and parameterizations may produce likelihood values on substantially different scales, the resulting scores are not directly comparable across edges\. We therefore standardize each edge independently using the median and median absolute deviation \(MAD\) of its training scores:

med\(t,p\)\\displaystyle\\text\{med\}\_\{\(t,p\)\}=medianm∈𝐗𝚝𝚛0​\(−ℓ​\(b\(t,p\)∗;𝐮m\)\),\\displaystyle=\\underset\{m\\in\\mathbf\{X\}\_\{\\mathtt\{tr\}\}^\{0\}\}\{\\text\{median\}\}\\big\(\-\\ell\(b^\{\*\}\_\{\(t,p\)\};\\mathbf\{u\}\_\{m\}\)\\big\),\(11\)MAD\(t,p\)\\displaystyle\\text\{MAD\}\_\{\(t,p\)\}=medianm∈𝐗𝚝𝚛0​\|−ℓ​\(b\(t,p\)∗;𝐮m\)−med\(t,p\)\|\.\\displaystyle=\\underset\{m\\in\\mathbf\{X\}\_\{\\mathtt\{tr\}\}^\{0\}\}\{\\text\{median\}\}\\big\|\{\-\\ell\(b^\{\*\}\_\{\(t,p\)\};\\mathbf\{u\}\_\{m\}\)\}\-\\text\{med\}\_\{\(t,p\)\}\\big\|\.\(12\)The median and MAD are preferred over the mean and standard deviation because of their robustness to outliers\. For a generic observation𝐮m\\mathbf\{u\}\_\{m\}, the standardized score at edge\(t,p\)\(t,p\)is defined as

𝚜𝚌𝚘𝚛𝚎​\(\(t,p\),𝐮m\)=−ℓ​\(b\(t,p\)∗;𝐮m\)−med\(t,p\)MAD\(t,p\),\\mathtt\{score\}\\big\(\(t,p\),\\mathbf\{u\}\_\{m\}\\big\)=\\frac\{\-\\ell\(b^\{\*\}\_\{\(t,p\)\};\\mathbf\{u\}\_\{m\}\)\-\\text\{med\}\_\{\(t,p\)\}\}\{\\text\{MAD\}\_\{\(t,p\)\}\},\(13\)which can be interpreted as the number of MAD units by which the edge score deviates from the ordinary baseline\. Collecting these values over all edges yields the score vector

𝐬m=\{𝚜𝚌𝚘𝚛𝚎​\(\(t,p\),𝐮m\):\(t,p\)∈ℰ→d−1\},\\mathbf\{s\}\_\{m\}=\\Big\\\{\\mathtt\{score\}\\big\(\(t,p\),\\mathbf\{u\}\_\{m\}\\big\):\(t,p\)\\in\\mathcal\{E\}\_\{\\to d\-1\}\\Big\\\},which summarizes the contribution of each pair\-copula to the anomaly assessment and forms the basis for both global detection and local interpretation\.

##### Global anomaly score\.

While various anomaly detection strategies can be applied to𝐬m\\mathbf\{s\}\_\{m\}, we employ a simple aggregation rule and define a scalar anomaly scoreSmS\_\{m\}as the average of theκ\\kappalargest edge values in𝐬m\\mathbf\{s\}\_\{m\}:

Sm=1κ​∑\(t,p\)∈𝒯κ\(m\)𝚜𝚌𝚘𝚛𝚎​\(\(t,p\),𝐮m\),S\_\{m\}=\\frac\{1\}\{\\kappa\}\\sum\_\{\(t,p\)\\in\\mathcal\{T\}\_\{\\kappa\}^\{\(m\)\}\}\\mathtt\{score\}\\big\(\(t,p\),\\mathbf\{u\}\_\{m\}\\big\),\(14\)where𝒯κ\(m\)\\mathcal\{T\}\_\{\\kappa\}^\{\(m\)\}denotes the set of theκ\\kappahighest\-scoring edges for observationmmandκ∈\{1,…,d​\(d−1\)/2\}\\kappa\\in\\\{1,\\dots,d\(d\-1\)/2\\\}is a user\-configurable parameter\.

##### Conformal anomaly detection\.

The aggregate scoreSmS\_\{m\}provides a scalar measure of abnormality\. We next convert this score into a statistically calibrated anomaly detector using split conformal prediction \(CP\)\[angelopoulos2023conformal\], introduced in Sect\.[2\.2](https://arxiv.org/html/2607.25020#S2.SS2)\.

We formulate anomaly detection as a binary classification problem, assigning label0to ordinary observations and label11to anomalous observations\. Given a test observation𝐱m\\mathbf\{x\}\_\{m\}with unknown labelymy\_\{m\}, letSmS\_\{m\}denote its aggregate score computed from the corresponding pseudo\-observation𝐮m\\mathbf\{u\}\_\{m\}\. We assess the compatibility ofSmS\_\{m\}with each class through class\-conditional non\-conformity functions

𝚗𝚌𝚏​\(Sm;0\)=Sm−𝚍𝚂,𝚗𝚌𝚏​\(Sm;1\)=𝚍𝚂−Sm,\\mathtt\{ncf\}\(S\_\{m\};0\)=S\_\{m\}\-\\mathtt\{d\_\{S\}\},\\qquad\\mathtt\{ncf\}\(S\_\{m\};1\)=\\mathtt\{d\_\{S\}\}\-S\_\{m\},where𝚍𝚂\\mathtt\{d\_\{S\}\}is a threshold separating ordinary and anomalous scores\. In our experiments,𝚍𝚂\\mathtt\{d\_\{S\}\}is chosen as the 95th percentile of the aggregate scores computed over𝐗𝚝𝚛0\\mathbf\{X\}\_\{\\mathtt\{tr\}\}^\{0\}\. The non\-conformity score is small when the observation is consistent with the corresponding class and large otherwise\.

Given a desired confidence levelα\\alpha, we adopt the Mondrian conformal framework\[vovk2005algorithmic\]and calibrate the two classes separately by applying𝚗𝚌𝚏\\mathtt\{ncf\}to the labeled observations in𝐗𝚌𝚊0\\mathbf\{X\}\_\{\\mathtt\{ca\}\}^\{0\}and𝐗𝚌𝚊1\\mathbf\{X\}\_\{\\mathtt\{ca\}\}^\{1\}, respectively\. The thresholdsτα0\\tau\_\{\\alpha\}^\{0\}andτα1\\tau\_\{\\alpha\}^\{1\}are then defined as the⌈\(n0\+1\)​\(1−α\)⌉/n0\\lceil\(n\_\{0\}\+1\)\(1\-\\alpha\)\\rceil/n\_\{0\}and⌈\(n1\+1\)​\(1−α\)⌉/n1\\lceil\(n\_\{1\}\+1\)\(1\-\\alpha\)\\rceil/n\_\{1\}empirical quantiles of the corresponding non\-conformity scores, wheren0=\|𝐗𝚌𝚊0\|n\_\{0\}=\|\\mathbf\{X\}\_\{\\mathtt\{ca\}\}^\{0\}\|andn1=\|𝐗𝚌𝚊1\|n\_\{1\}=\|\\mathbf\{X\}\_\{\\mathtt\{ca\}\}^\{1\}\|\. The*prediction region*associated with𝐱m\\mathbf\{x\}\_\{m\}is:

Γα​\(𝐱m\)=\{y∈\{0,1\}:𝚗𝚌𝚏​\(Sm,y\)≤ταy\}\.\\Gamma\_\{\\alpha\}\(\\mathbf\{x\}\_\{m\}\)=\\big\\\{y\\in\\\{0,1\\\}:\\mathtt\{ncf\}\(S\_\{m\},y\)\\leq\\tau\_\{\\alpha\}^\{y\}\\big\\\}\.\(15\)An observation receives label0when𝚗𝚌𝚏​\(Sm,0\)≤τα0\\mathtt\{ncf\}\(S\_\{m\},0\)\\leq\\tau\_\{\\alpha\}^\{0\}and label11when𝚗𝚌𝚏​\(Sm,1\)≤τα1\\mathtt\{ncf\}\(S\_\{m\},1\)\\leq\\tau\_\{\\alpha\}^\{1\}\. Consequently, the prediction region can be\{0\}\\\{0\\\},\{1\}\\\{1\\\},\{0,1\}\\\{0,1\\\}, or∅\\emptyset, corresponding respectively to a unique ordinary prediction, a unique anomalous prediction, ambiguity between the two classes, or incompatibility with both\.

Under the assumption of exchangeability between calibration and test observations belonging to the same class, Mondrian conformal prediction*guarantees class\-conditional coverage*:

ℙ\(𝐱m,ym\)​\(ym∈Γα​\(𝐱m\)\)≥1−α,\\mathbb\{P\}\_\{\(\\mathbf\{x\}\_\{m\},y\_\{m\}\)\}\\bigl\(y\_\{m\}\\in\\Gamma\_\{\\alpha\}\(\\mathbf\{x\}\_\{m\}\)\\bigr\)\\geq 1\-\\alpha,meaning that the prediction region contains the true label with probability at least1−α1\-\\alphafor each class\[vovk2005algorithmic\]\. This distinction is particularly important in anomaly detection, where the anomalous class is often heavily underrepresented\. Standard conformal prediction guarantees only marginal coverage over the entire population, which may be dominated by ordinary observations and therefore provides substantially weaker guarantees for anomalies\. By calibrating each class separately, Mondrian conformal prediction ensures the desired coverage level for both ordinary and anomalous observations\.

##### Interpretability\.

The standardized edge scores naturally provide an interpretable explanation of detected anomalies\. Each pair\-copulaCi​j∣DC\_\{ij\\mid D\}models the dependence between variablesiiandjjconditional on a setDD\. Consequently, a large score at the associated edge indicates that the anomaly is associated with an unexpected dependence pattern in that specific relationship\. Rather than merely detecting the presence of an anomaly, the proposed framework identifies which variable interactions contribute most strongly to the anomalous behavior\.

## 4Experiments

##### Pair\-copula implementations\.

To support gradient\-based fitting, we developed a fully differentiable PyTorch implementation of the Gaussian and Student\-t elliptical copulas, together with the Clayton, Frank, Gumbel, and Joe Archimedean families\. Copula densities and h\-functions are implemented in a differentiable and numerically stable manner, with real\-valued parameters mapped to their valid domains through scaled sigmoid transformations\. The supported families span a broad range of dependence structures; parameter bounds, tail\-dependence properties and closed\-form expressions with illustrative examples are provided in AppendixLABEL:app:families\. Numerical stability is further ensured through log\-space density evaluation, input clamping away from 0 and 1, and masking of non\-finite likelihood contributions\.

##### Experimental setup\.

We evaluate the proposed methodology on two datasets: the Wilt dataset from the ADBench benchmark\[han2022adbench\]and a proprietary real\-world dataset of operational measurements from a sewerage network in northern Italy\. Both datasets are well suited to copula\-based modeling, as they consist exclusively of continuous variables, contain a sufficient number of observations and labeled anomalies to support both training and evaluation, and involve a relatively small number of features, resulting in compact D\-vine decompositions\.

Unless otherwise stated, we report results for the penalized fitting procedure introduced in Eq\.[10](https://arxiv.org/html/2607.25020#S3.E10)and the corresponding anomaly detection framework\. Comparisons with the non\-penalized variant are discussed where relevant, while complete results are provided in AppendixLABEL:app:results\.

For both datasets, the candidate family set comprises the Gaussian, Student\-t, Clayton, Joe, Frank, and Gumbel copulas\. During fitting, each pair\-copula is optimized for 250 epochs in the initial edge\-wise stage, followed by 200 epochs of joint refinement for each surviving beam state inℬd−1\\mathcal\{B\}\_\{d\-1\}\. The Clarke test is used to evaluate competing edge configurations with the maximum branching factor set toωmax=4\\omega\_\{\\max\}=4, while the beam width is limited toβmax=8\\beta\_\{\\max\}=8\. Unless otherwise stated, the anomaly penalization coefficient in \([10](https://arxiv.org/html/2607.25020#S3.E10)\) is fixed atλ=0\.1\\lambda=0\.1\.

![Refer to caption](https://arxiv.org/html/2607.25020v1/Figures/44_WILT_NEW_OPO.png)Figure 3:Wilt:Pairwise scatterplots of the ordered pseudo\-observations\. Variables pairs plotted together form the inputs to the pair\-copulas at the first treeT1T\_\{1\}\.![Refer to caption](https://arxiv.org/html/2607.25020v1/Figures/44_WILT_NEW_FittedCopula.png)Figure 4:Wilt:Fitted pair\-copulas of the D\-vine, shown as joint densities\. Each panel can be read in isolation, as the dependence accumulated in earlier trees is embedded in the inputs of each panel through the h\-function cascade\.
### 4\.1Wilt Dataset Results

The Wilt dataset is derived from a remote\-sensing task in which satellite image segments of forested areas are classified as either diseased trees or other land cover, with the diseased\-tree segments treated as anomalies\. It is comprised of roughly 4500 ordinary and 250 anomalous observation across five spectral features\. After applying the estimated marginal CDFs and variable ordering, the resulting pseudo\-observations reveal two different dependence patterns for ordinary and anomalous observations, as can be seen in the third panel of Fig\.[3](https://arxiv.org/html/2607.25020#S4.F3), making the dataset an ideal candidate for our copula\-based approach\.

Fig\.[4](https://arxiv.org/html/2607.25020#S4.F4)shows the D\-vine configuration estimated by the proposed fitting procedure\. Each panel visualizes the joint density induced by the fitted pair\-copula combined with standard Gaussian marginals\. Since the true marginal distributions of conditional variables are generally unavailable beyond the first tree, Gaussian marginals are adopted solely for visualization purposes\. This choice does not affect the interpretation of the dependence structure, as copulas are invariant to the marginal distributions\[nelsen2006introduction\]\.

We notice that the penalization does not affect all edges uniformly as most retain similar family parameter configurations as in the non penalized fit \(AppendixLABEL:app:results\), with its effects being concentrated on edges\(1,3\)\(1,3\)and\(2,2\)\(2,2\)\. Under penalization, edge\(1,3\)\(1,3\)switches from a Gumbel to a Joe copula, with the latter exhibiting stronger upper tail dependence\. Edge\(2,2\)\(2,2\)maintains a Student\-t copula but the value of the second parameter, the degrees of freedom, drops substantially, indicating heavier tails and a departure from Gaussian behavior\. As can be seen in the third panel of Fig\.[3](https://arxiv.org/html/2607.25020#S4.F3)and in Fig\.[4](https://arxiv.org/html/2607.25020#S4.F4),\(1,3\)\(1,3\)is the edge inT1T\_\{1\}with the clearest separation between ordinary and anomalous observations and, as the inputs of\(2,2\)\(2,2\)are derived in part from\(1,3\)\(1,3\)through the h\-function cascade, this separation propagates to the second tree\.

![Refer to caption](https://arxiv.org/html/2607.25020v1/Figures/44_WILT_NEW_EDGESCORES_with_pen_attribution.png)Figure 5:Wilt:Standardized edge scores𝚜𝚌𝚘𝚛𝚎​\(\(t,p\),𝐮m\)\\mathtt\{score\}\\big\(\(t,p\),\\mathbf\{u\}\_\{m\}\\big\)for the Wilt dataset across the calibration and test sets, for fitting with penalization\.The relevance of these edges in identifying potential anomalies is also visible in Fig\.[5](https://arxiv.org/html/2607.25020#S4.F5), which shows the standardized scores𝐬m\\mathbf\{s\}\_\{m\}for each edge in the decomposition and for each observationmmin𝐗𝚌𝚊0\\mathbf\{X\}\_\{\\mathtt\{ca\}\}^\{0\},𝐗𝚌𝚊1\\mathbf\{X\}\_\{\\mathtt\{ca\}\}^\{1\},𝐗𝚝𝚎0\\mathbf\{X\}\_\{\\mathtt\{te\}\}^\{0\}and𝐗𝚝𝚎1\\mathbf\{X\}\_\{\\mathtt\{te\}\}^\{1\}\. Edges\(1,3\)\(1,3\)and\(2,2\)\(2,2\)present significantly higher scores for observations belonging to the anomalous sets𝐗𝚌𝚊1\\mathbf\{X\}\_\{\\mathtt\{ca\}\}^\{1\}and𝐗𝚝𝚎1\\mathbf\{X\}\_\{\\mathtt\{te\}\}^\{1\}compared to those belonging to their ordinary counterparts, confirming that these edges model the variable relationships through which anomalies can be detected\. Lastly, it can be noted that the calibration and test splits exhibit consistent patterns within each class, indicating that the standardized edge scores generalize from the calibration set to unseen observations\. For this dataset, we compute the aggregate anomaly score following Eq\.[14](https://arxiv.org/html/2607.25020#S3.E14)withκ=2\\kappa=2\. After setting the significance levelα=0\.1\\alpha=0\.1, we calibrate the two class\-conditional thresholdsτα0\\tau^\{0\}\_\{\\alpha\}andτα1\\tau^\{1\}\_\{\\alpha\}on𝐗𝚌𝚊0\\mathbf\{X\}\_\{\\mathtt\{ca\}\}^\{0\}and𝐗𝚌𝚊1\\mathbf\{X\}\_\{\\mathtt\{ca\}\}^\{1\}as per Sect\.[3\.2](https://arxiv.org/html/2607.25020#S3.SS2)and assign each observation in𝐗𝚝𝚎0\\mathbf\{X\}\_\{\\mathtt\{te\}\}^\{0\}and𝐗𝚝𝚎1\\mathbf\{X\}\_\{\\mathtt\{te\}\}^\{1\}a predictive regionΓα​\(𝐱m\)\\Gamma\_\{\\alpha\}\(\\mathbf\{x\}\_\{m\}\)\. Table[1](https://arxiv.org/html/2607.25020#S4.T1)reports the class conditional conformal prediction results on the Wilt dataset forα=0\.1\\alpha=0\.1, with both the ordinary and anomalous classes reaching the target coverage level of0\.90\.9\. Most observations receive a single correct label, none are assigned an empty set, and only a minority receive an ambiguous prediction\. Fig\.[6](https://arxiv.org/html/2607.25020#S4.F6)provides a graphical view of the same result over the two most discriminative edges,\(1,3\)\(1,3\)and\(2,2\)\(2,2\)\. Correct single\-label assignments concentrate where the two classes are well separated, while ambiguous or incorrect assignments tend to occur in the intermediate region where ordinary and anomalous observations overlap\.

![Refer to caption](https://arxiv.org/html/2607.25020v1/Figures/WILT_SelectedEdges.png)Figure 6:Wilt:Anomaly detection results on the most relevant edges; color indicates prediction correctness, shape distinguishes ordinary from anomalous observations\.![Refer to caption](https://arxiv.org/html/2607.25020v1/Figures/BRIANZA_NEW_EDGESCORES_with_pen_attribution.png)Figure 7:Sewerage:Standardized edge scores𝚜𝚌𝚘𝚛𝚎​\(\(t,p\),𝐮m\)\\mathtt\{score\}\\big\(\(t,p\),\\mathbf\{u\}\_\{m\}\\big\)for the Sewerage dataset across the calibration and test sets, for fitting with penalization\.
### 4\.2Sewerage Dataset Results

Our second dataset is a multivariate time\-series of real\-world operational measurements from a sewerage network in northern Italy, containing four months of hourly level recordings from four sensors placed at neighbouring points within the same network, referred to as Sensors 1\-4\. Each recording is treated as an individual observation\. Under ordinary dry\-weather conditions, all four sensors follow the same regular daily pattern, and rainfall events produce simultaneous spikes in level across all four\. During the observation period, Sensor 4 enters an anomalous state in which its readings, while remaining within an ordinary operational range, no longer follow the behavior shared by the rest of the network\.

Following the same procedure as for the Wilt dataset, we fit a D\-vine for the four sensors and report the resulting standardized edge scores𝐬m\\mathbf\{s\}\_\{m\}in Fig\.[7](https://arxiv.org/html/2607.25020#S4.F7)\. Edges\(1,2\)\(1,2\),\(1,3\)\(1,3\)and\(3,1\)\(3,1\)clearly discriminate between ordinary and anomalous observations, with significantly higher values on the anomalous sets\. Notably, after accounting for variable ordering,\(1,2\)\(1,2\)and\(1,3\)\(1,3\)are the pair\-copulas inT1T\_\{1\}taking as input Sensor 4, demonstrating the method’s ability to identify the cause of the deviation from expected behavior\. Lastly, we compute the aggregate score of Eq\.[14](https://arxiv.org/html/2607.25020#S3.E14)withκ=3\\kappa=3and apply the same conformal procedure atα=0\.1\\alpha=0\.1, obtaining the predictive regionsΓα\\Gamma\_\{\\alpha\}for all observations in𝐗𝚝𝚎0\\mathbf\{X\}\_\{\\mathtt\{te\}\}^\{0\}and𝐗𝚝𝚎1\\mathbf\{X\}\_\{\\mathtt\{te\}\}^\{1\}whose results are reported in Table[1](https://arxiv.org/html/2607.25020#S4.T1)and Fig\.[8](https://arxiv.org/html/2607.25020#S4.F8)\. Both classes reach the desired coverage level and, as for the Wilt dataset, few observations are misclassified or assigned an ambiguous prediction, with none receiving an empty set\. Correct single label assignment concentrate where the two classes are clearly separated\.

Table 1:Class\-conditional conformal prediction results on both datasets atα=0\.1\\alpha=0\.1\. Coverage is the fraction of observations whose prediction region contains the true label\.![Refer to caption](https://arxiv.org/html/2607.25020v1/Figures/BRIANZA_SelectedEdges.png)Figure 8:Sewerage:Anomaly detection results on the most relevant edges; color indicates prediction correctness, shape distinguishes ordinary from anomalous observations\.

## 5Conclusion

We presented a novel framework for localized anomaly detection based on differentiable D\-vine copulas\. At its core is a beam\-search fitting procedure that maintains multiple candidate D\-vine configurations throughout the construction process, avoiding the myopic decisions of standard greedy methods\. This is enabled by a fully differentiable PyTorch implementation of several pair\-copula families, allowing gradient\-based estimation and joint refinement of full D\-vine configurations\. Building on the fitted D\-vine, we introduced an anomaly detection approach that exploits the hierarchical pair\-copula decomposition to produce global anomaly scores and localized explanations\. Class\-conditional statistical guarantees are provided via Mondrian conformal prediction, enabling reliable uncertainty quantification even in unbalanced anomaly detection settings\.

Experiments on both a public benchmark dataset and a real\-world sewerage monitoring application demonstrate the effectiveness of the proposed framework\. Beyond accurate anomaly detection, the standardized edge scores consistently identify the variable relationships responsible for anomalous behavior, providing interpretable insights into the underlying causes of detected deviations\. Future work will explore alternative vine structures, such as C\-vines, which may provide additional insight into hierarchical and potentially causal relationships among variables by explicitly modeling the influence of central variables on remaining system components\. We also plan to extend the framework to temporal and streaming settings, where dependence patterns evolve over time\.

\{credits\}

#### 5\.0\.1\\discintname

The authors have no competing interests to declare that are relevant to the content of this article\.

## References