# Quasi-Bayes meets Vines

David Huk
Department of Statistics
University of Warwick
<EMAIL_ADDRESS>
&Yuanhe Zhang
Department of Statistics
University of Warwick
<EMAIL_ADDRESS>
Mark Steel
Department of Statistics
University of Warwick
<EMAIL_ADDRESS>
&Ritabrata Dutta
Department of Statistics
University of Warwick
<EMAIL_ADDRESS>

###### Abstract

Recently proposed quasi-Bayesian (QB) methods (Fong et al., 2023) initiated a
new era in Bayesian computation by directly constructing the Bayesian
predictive distribution through recursion, removing the need for expensive
computations involved in sampling the Bayesian posterior distribution. This
has proved to be data-efficient for univariate predictions, but extensions to
multiple dimensions rely on a conditional decomposition resulting from
predefined assumptions on the kernel of the Dirichlet Process Mixture Model,
which is the implicit nonparametric model used. Here, we propose a different
way to extend Quasi-Bayesian prediction to high dimensions through the use of
Sklar’s theorem by decomposing the predictive distribution into one-
dimensional predictive marginals and a high-dimensional copula. Thus, we use
the efficient recursive QB construction for the one-dimensional marginals and
model the dependence using highly expressive vine copulas. Further, we tune
hyperparameters using robust divergences (eg. energy score) and show that our
proposed Quasi-Bayesian Vine (QB-Vine) is a fully non-parametric density
estimator with _an analytical form_ and convergence rate independent of the
dimension of data in some situations. Our experiments illustrate that the QB-
Vine is appropriate for high dimensional distributions ($\sim$64), needs very
few samples to train ($\sim$200) and outperforms state-of-the-art methods with
analytical forms for density estimation and supervised tasks by a considerable
margin.

##  Introduction

The estimation of joint densities $\mathbf{p}(\mathbf{x})$ is a cornerstone of
machine learning as a looking glass into the underlying data-generating
process of multivariate data $\mathbf{x}$. Methods that support explicit
density evaluation are crucial in probabilistic modelling, with applications
in variational methods (Kingma and Welling, 2013; Rezende et al., 2014),
Importance Sampling (Arbel et al., 2021; Matthews et al., 2022), Sequential
Monte Carlo (Gu et al., 2015), Markov Chain Monte Carlo (MCMC) (Song et al.,
2017; Hoffman et al., 2019) and simulation-based inference (Papamakarios et
al., 2019; Lueckmann et al., 2021). A prominent example are Normalising Flows
(NF) Papamakarios et al. (2017); Durkan et al. (2019); Papamakarios et al.
(2021), leveraging deep networks with invertible transformations for
analytical expressions and sampling. Despite impressive performances, they
require meticulous manual hyperparameter tuning and large amounts of data to
train. Bayesian methods are another attractive approach for analytical density
modelling where the central object of interest is the predictive density
$\mathbf{p}^{(n)}(\mathbf{x})$, with the Dirichlet Process Mixture Model
(DPMM) Hjort et al. (2010) as the canonical nonparametric choice. However,
obtaining analytical predictive densities commonly relies on computationally
expensive MCMC methods scaling poorly to high-dimensions111For reference, with
$d$ the dimensionality of the parameter space, Random Walk Metropolis-Hastings
scales like $\mathcal{O}(d^{2})$ (Roberts et al., 1997), Metropolis-adjusted
Langevin algorithm like $\mathcal{O}(d^{5/4})$ (Roberts and Rosenthal, 1998),
and Hamiltonian Monte Carlo like $\mathcal{O}(d^{4/3})$ (Beskos et al., 2013).

Recently, Fong et al. (2023) heralded a revolution in Bayesian methods with
efficient constructions of Quasi-Bayesian (QB) sequences of predictive
densities through recursion alone, liberating Bayes from the reliance on MCMC.
The seminal univariate Recursive Bayesian Predictive (R-BP) (Hahn et al.,
2018), its multivariate extension (Fong et al., 2023) and the AutoRegressive
Bayesian Predictive (AR-BP) (Ghalebikesabi et al., 2022) are all QB models
targeting the predictive mean of the DPMM, with an analytical recursive form
driven by bivariate copula updates. Notably, the AR-BP demonstrated superior
density estimation capabilities compared to a suite of competitors on varied
supervised and unsupervised datasets, is orders of magnitude faster than
standard Bayesian methods and is data-efficient, i.e. does not require large
amounts of training data to be effective. In multivariate QB predictives,
analytical expressions are derived by enforcing assumptions on the DPMM kernel
structure. The kernel of the R-BP is set to be independent across dimensions,
where this strong assumption is too constraining for more complex data. This
is relaxed in the AR-BP through an autoregressive form of the kernel, where
each kernel mean is a similarity function of previous dimensions modelled
through kernels or deep autoregressive networks, depending on the complexity
of the data. The former relies on a fixed form that is often too simplistic
while the latter loses the appeal of a data-efficient predictive like in the
R-BP.

In this paper, we introduce the Quasi-Bayesian Vine (QB-Vine) for density
estimation and supervised learning, obtained by applying Sklar’s theorem
(Sklar, 1959) to the joint predictive, thereby dissecting it into univariate
marginal predictive densities and a multivariate copula. Marginal predictives
are modelled with the data-efficient univariate R-BP while the multivariate
copula is modelled with a highly flexible vine copula, see the diagram in
Figure 1.

Figure 1: Overview of the Quasi-Bayesian Vine model. The joint predictive
density $\mathbf{p}^{(n)}\left(x_{1},\ldots,x_{d}\right)$ is modelled
nonparametrically through a copula decomposition into univariate marginal
predictives and a multivariate copula. Left: Univariate predictive densities
$p_{i}^{(n)}$ are modelled with a data-efficient Quasi-Bayesian recursion
using observed samples $\\{x_{i}^{(k)}\\}_{k=1}^{n}$ alone, bypassing
posterior integration steps. Right: A vine copula model decomposes the high-
dimensional copula into $\frac{d(d-1)}{2}$ two-dimensional copulas
nonparametrically estimated with kernel density copula estimators represented
by a graphical structure between dimensions

The main contributions of our work are as follows:

  * •

This decomposition frees us from the need to assume specific kernel structures
of the DPMM, but preserves the data-efficiency of QB methods while making it
more effective on high-dimensional data.

  * •

Under the assumption of well-identified simplified vine copula model, we show
that the QB-Vine attains a convergence rate that is independent of dimension.

  * •

The decomposition of the joint density and use of energy score to tune
hyperparameters makes our algorithm amenable to efficient parallelisation with
significant computational gains.

Our paper is structured as follows. In Section 2 we introduce Quasi-Bayesian
prediction, recapitulating the R-BP construction. In Section 3 we formulate
the Quasi-Bayesian Vine model. We provide a succinct survey of related work in
Section 4 and compare related methods to the QB-Vine in Section 5 on a range
of datasets for density estimation, regression and classification, achieving
state-of-the-art performance with our model. We conclude with a discussion in
Section 6.

##  Quasi-Bayesian prediction

#### Notation.

Let $p(\mathbf{x})$ be a multivariate probability density function over
$\mathcal{X}\subseteq\mathbb{R}^{d}$, from which we observe i.i.d. samples
$\mathcal{D}_{P}=\\{\mathbf{x}^{k}\\}_{k=1}^{K}\sim p(\mathbf{x})$. Similarly,
let $p_{1}(x_{1}),\ldots,p_{d}(x_{d})$ be the marginal densities of
$p(\mathbf{x})$, each over (a subset of) $\mathbb{R}$ with corresponding
i.i.d. samples $\mathcal{D}_{P_{i}}=\\{{x_{i}}^{k}\\}_{k=1}^{K}\sim
p_{i}({x_{i}})$ for $i=1,\ldots,d$, and assumed to all be continuous. Further,
let $P$ and $P_{1},\ldots,P_{d}$ be the respective cumulative distribution
functions (cdf’s) of the previously mentioned densities, in order of
appearance. Finally, when discussing predictive densities, we will use a
superscript, plain for data (e.g. $x^{n}$) and in parentheses for functions
(e.g. $p^{(n)}$) to indicate the predictive at the $n^{\text{th}}$ step,
distinguishing them from the subscripts kept for dimension.

#### Bayesian predictive densities as copula updates.

Consider the univariate Bayesian predictive density $p^{(n)}$ for a future
observation $x$ given seen i.i.d. data $x^{1:n}\in\mathbb{R}$ with a
likelihood $f$ of the data and a posterior $\pi^{(n)}$ for the model
parameters $\theta$ after $n$ observations. By Bayes rule, we have:

$p^{(n)}(x|x^{1:n})=\frac{\int f(x|\theta)\cdot
f(x^{n}|\theta)\cdot\pi^{(n-1)}(\theta|x^{1:n-1})~{}d\theta}{p^{(n-1)}(x^{n}|x^{1:n-1})}.$
(1)

As identified by Hahn et al. (2018); Fong et al. (2023); Holmes and Walker
(2023), multiplying and dividing by the predictive from the previous step
$p^{(n-1)}$, we arrive at

$\begin{split}p^{(n)}(x|x^{1:n})&=p^{(n-1)}(x|x^{1:n-1})\cdot\frac{\overbrace{\int
f(x|\theta)\cdot
f(x^{n}|\theta)\cdot\pi^{(n-1)}(\theta|x^{1:n-1})~{}d\theta}^{\text{ Joint
density for $x,x^{n}$}}}{\underbrace{p^{(n-1)}(x^{n}|x^{1:n-1})}_{\text{
Marginal for $x^{n}$}}\cdot\underbrace{p^{(n-1)}(x|x^{1:n-1})}_{\text{
Marginal for $x$}}}\\\ &=p^{(n-1)}(x|x^{1:n-1})\cdot
c^{(n)}\left(P^{(n-1)}(x),P^{(n-1)}(x^{n})\right)\end{split}$ (2)

which identifies the involvement of copulas (see Appendix A, Elidan (2013);
Größer and Okhrin (2022)) in this recursive equation for the predictive
densities. The second term on the right-hand side of (2) is seen to be a
symmetric bivariate copula density function with the property that
$c^{(n)}(\cdot,\cdot)\rightarrow 1$ to ensure $p^{(n)}$ converges
asymptotically with $n$. It can be shown that every univariate Bayesian model
can be written in this form Hahn et al. (2018), and so has a unique copula
sequence characterising its predictive updates by de Finetti’s theorem (De
Finetti, 1937; Hewitt and Savage, 1955). The choice of the predictive sequence
then corresponds to an implicit choice of likelihood and prior (Berti et al.,
2023; Garelli et al., 2024).

#### Recursive Bayesian Predictives.

Due to the integration over the posterior density often being intractable in
practice, identifying the update copulas $c^{(n)}$ analytically is generally
impossible. Therefore, Hahn et al. (2018) propose a nonparametric density
estimator termed recursive Bayesian predictive (R-BP) as a Dirichlet Process
Mixture Model (DPMM) inspired recursion emulating (2). They derive the correct
Bayesian update under a DPMM for step $1$ and use it for all future steps
$m>1$ (derivations shown in Appendix B.2). The update copula of the R-BP is a
mixture between the independent and Gaussian copula, thereby deviating from
the true (unknown) form of the Bayesian recursion copulas for a DPMM. For an
initial choice of predictive density $p^{(0)}$ and distribution $P^{(0)}$, the
obtained analytical expression for the R-BP recursion has the following
predictive density

$p^{(n)}(x)=p^{(n-1)}(x)\cdot\left[(1-\alpha_{n})+\alpha_{n}\cdot
c_{\rho}(P^{(n-1)}(x),P^{(n-1)}(x^{n}))\right]$ (3)

with $c_{\rho}\in[0,1]$ being a bivariate Gaussian copula with covariance
$\rho$. The corresponding cdf also admits an analytical expression as follows

$P^{(n)}(x)=(1-\alpha_{n})\cdot
P^{(n-1)}(x)+\alpha_{n}H_{\rho}(P^{(n-1)}(x),P^{(n-1)}(x^{n})).$ (4)

Here, for $\Phi$ the standard univariate Gaussian distribution,

$H_{\rho}(u,v)=\Phi\left(\frac{\Phi^{-1}(u)-\rho\cdot\Phi^{-1}(v)}{\sqrt{1-\rho^{2}}}\right)$
(5)

is a conditional Gaussian copula distribution with covariance $\rho$ treated
as a hyperparameter and where $\alpha_{k}=(2-\frac{1}{k})\frac{1}{k+1}$ is a
sequence of weights converging to $0$ (See supplement E of Fong et al. (2023)
for a more detailed explanation of the weights). This univariate R-BP model
has been shown to converge to a limiting distribution $P^{(\infty)}$ with
density $p^{(\infty)}$ termed the _martingale posterior_ (Fong et al., 2023;
Fortini and Petrone, 2023; Holmes and Walker, 2023; Berti et al., 2006;
Garelli et al., 2024) which is a quasi-Bayesian object describing current
uncertainty (See Appendix B for an introduction). Frequentist results such as
Kullback-Leibler convergence to the true density are given in Hahn et al.
(2018). An extension to the multivariate case was studied by Fong et al.
(2023) and refined by Ghalebikesabi et al. (2022) for data with more complex
dependence structures. The recursive formulation of the R-BP with an
analytical form ensures fast updates of predictives whilst the use of copulas
bypasses the need to evaluate a normalising constant (as is the case in
Newton’s algorithm (Newton et al., 1998; Newton, 2002)). Consequently, this
approach is free from the reliance on MCMC to approximate a posterior density,
making it much faster than regular DPMMs. While this formulation does not
correspond to a Bayesian model, as argued by Berti et al. (2004, 2023); Fong
et al. (2023); Holmes and Walker (2023), if the recursive updates are
conditionally identically distributed (as is the case for the recursion of
(4)), they still exhibit desirable Bayesian characteristics such as coherence,
regularization, and asymptotic exchangeability, motivating the Quasi-Bayesian
name as used in Fortini and Petrone (2020).

##  Quasi-Bayesian Vine prediction

We propose the Quasi-Bayesian Vine (QB-Vine) for efficiently modeling a high-
dimensional predictive density $\mathbf{p}^{(n)}$ (and distribution
$\mathbf{P}^{(n)}$). The efficiency is achieved by splitting the joint density
into predictive marginals for each dimension and a high-dimensional copula;
then using two different nonparametric density estimation schemes for the
marginals and the vine, correspondingly the Quasi-Bayes recursion (described
in Equation 3 and 4) for each of the marginals and a vine copula to model the
high-dimensional dependency. We first start by describing how we decompose the
joint predictive density adapting Sklar’s theorem (Sklar (1959), see Appendix
A),

######  Theorem 3.1.

Let $\mathbf{P}^{(n)}$ be an $d$-dimensional predictive distribution function
with continuous marginal distributions
$P^{(n)}_{1},P^{(n)}_{2},\ldots,P^{(n)}_{d}$. Then there exists a copula
distribution $\mathbf{C^{(n)}}$ such that for all
$\mathbf{x}=(x_{1},x_{2},\ldots,x_{d})\in\mathbb{R}^{d}$:

$\mathbf{P}^{(n)}(x_{1},\ldots,x_{d})=\mathbf{C}^{(n)}(P_{1}^{(n)}(x_{1}),\ldots,P_{d}^{(n)}(x_{d}))$
(6)

And if a probability density function (pdf) is available:

$\mathbf{p}^{(n)}(x_{1},\ldots,x_{d})=p_{1}^{(n)}(x_{1})\cdot\ldots\cdot
p_{d}^{(n)}(x_{d})\cdot\mathbf{c}^{(n)}(P_{1}^{(n)}(x_{1}),\ldots,P_{d}^{(n)}(x_{d}))$
(7)

where $p_{1}^{(n)}(x_{1}),\ldots,p_{d}^{(n)}(x_{d})$ are the marginal
predictive probability density functions (pdf), and
$\mathbf{c}^{(n)}:[0,1]^{d}\rightarrow\mathbb{R}$ is the copula pdf.

The decomposition of (7) has two consequences for modelling
$\mathbf{p}^{(n)}$. Firstly, as the predictives $p_{i}^{(n)}$ are
unconditional marginal densities, their estimation can be done independently.
Secondly, by applying (7) to two consecutive predictive densities
$\mathbf{p}^{(m-1)}$ and $\mathbf{p}^{(m)}$, we obtain a recursive update for
predictive densities with two parts, of the form:

$\frac{\mathbf{p}^{(m)}}{\mathbf{p}^{(m-1)}}=\underbrace{\prod_{i=1}^{d}\Bigl{\\{}\frac{p_{i}^{(m)}}{p_{i}^{(m-1)}}\Bigr{\\}}}_{\text{Independent
recursions}}\cdot\underbrace{\frac{\mathbf{c}^{(m)}\Bigl{(}P_{1}^{(m)}(x_{1}),\ldots,P_{d}^{(m)}(x_{d})\Bigr{)}}{\mathbf{c}^{(m-1)}\Bigl{(}P_{1}^{(m-1)}(x_{1}),\ldots,P_{d}^{(m-1)}(x_{d})\Bigr{)}}}_{\text{Implicit
recursion on copulas}}\,.$ (8)

This decomposition fruitfully isolates updates to marginal predictive
densities from updates to their dependence structure, allowing us to model
each recursion separately; the marginal predictives follow a univariate
recursion _a la_ (2) while the copulas are free to follow another recursive
form. Further, as we are only interested in the joint predictive
$\mathbf{p}^{(n)}$, once marginal predictives are obtained, we only need to
fit a single copula $\mathbf{c}^{(n)}$ at the final iteration $n$ to recover
the joint predictive through (7). Crucially, we do not need an analytical form
for copula updates, leaving them as an implicit recursion on copula densities.
Motivated by this insight, we formulate the Quasi-Bayesian Vine as an
estimation problem with two parts, the recursive estimation of marginal
predictive densities and the estimation of a vine copula to model the high-
dimensional dependency.

#### Marginal predictive density estimation

We model marginal predictive densities ${p_{i}}^{(n)}(x_{i})$ and
distributions ${P_{i}}^{(n)}(x_{i})$, $\forall i=1,\ldots,d$ independently of
each other. We use the univariate R-BP approach described in Section 2 to
recursively obtain the analytical expression for both. For each dimension
separately, starting with an initial density $p^{(0)}_{i}$ and distribution
$P^{(0)}_{i}$, we follow the updates (3) for the density and (4) for the
distribution.

#### Simplified vine copulas for high-dimensional dependence

After estimating marginal predictives, we model the joint density of
$\left(u_{1}:={P_{1}}^{(n)}(x_{1}),\ldots,u_{d}:={P_{d}}^{(n)}(x_{d})\right)$
with a multivariate copula. We consider a highly flexible vine copula which
decomposes the joint copula density $c(u_{1},\ldots,u_{d})$ into
$\frac{d(d-1)}{2}$ bivariate conditional copulas to capture the dependence
structure.

Vine copulas are a class of copulas that provide a divide-and-conquer approach
to high-dimensional modelling by decomposing the joint copula density into
$\frac{d(d-1)}{2}$ bivariate copula terms. The main ingredient of a vine
copula decomposition is the following identity as a consequence of (7):

$p_{a|b}(x_{a}|x_{b})=c_{a,b}(P_{a}(x_{a}),P_{b}(x_{b}))\cdot p_{a}(x_{a})$
(9)

where $a,b$ are subsets of dimensions from $\\{1,\ldots,d\\}$. Vine copulas
rely on a conditional factorisation
$\mathbf{p}(x_{1},\ldots,x_{d})=\prod_{i=1}^{d}p_{i|<i}(x_{i}|x_{<i})$ to
which they repeatedly apply (9), rewriting any conditional densities as
copulas, thereby splitting the joint density into the $d$ marginal densities
and $\frac{d(d-1)}{2}$ bivariate copulas called pair copulas. The pair copulas
for each $i\neq j\in\\{,\ldots,d\\}$, take as input pairs of conditional
distributions
$(P_{i|\mathcal{S}_{ij}}(x_{i|\mathcal{S}_{ij}}),P_{j|\mathcal{S}_{ij}}(x_{j|\mathcal{S}_{ij}}))$
where
$\mathcal{S}_{ij}\subseteq\\{1,\ldots,d\\}\setminus\\{i,j\\}\cup\emptyset$ is
decided by the choice of the vine. A vine copula model thus has the form

${\mathbf{c}}(u_{1},\ldots,u_{d})=\prod_{i\neq
j}^{d(d-1)/2}{c}_{ij}(P_{i|\mathcal{S}_{ij}}(x_{i|\mathcal{S}_{ij}}),P_{j|\mathcal{S}_{ij}}(x_{j|\mathcal{S}_{ij}})|\mathcal{S}_{ij}).$
(10)

We note that these pair copulas start as unconditional bivariate copulas and
later capture higher orders of multivariate dependence by conditioning on the
set $\mathcal{S}$ itself. This decomposition is valid but only unique up to
the permutation of indexes. We provide an example of a three-dimensional vine
copula decomposition and an overview in Appendix A.3, referring the reader to
Czado (2019); Czado and Nagler (2022) for an introduction. In practice, we use
a _simplified vine copula_ model (Nagler and Czado, 2016; Nagler et al., 2017)
which removes the conditional dependence of each of the copula $c_{ij}$ on
$\mathcal{S}_{ij}$. This is an approximation reducing the complexity of the
model for dependency structure,

but provides a significant computation efficiency by reducing the size of tree
space and Nagler and Czado (2016) provides a dimension-independent convergence
rate when the simplified vine assumption is true, making simplified vine
copulas greatly appealing for high-dimensional models. Under this simplified
vine copula model assumption, we use nonparametric Kernel Density Estimator
(KDE) for each of the bivariate unconditional pair-copulas222For further
detail on KDE-based vine copulas, see Appendices A.2, A.3 and Nagler and Czado
(2016); Nagler et al. (2017).. Then ${c}_{ij}$ becomes a two-dimensional KDE
copula with the following expression:

${c}_{ij}(u,v)\,=\,\frac{\sum_{k=1}^{K}\phi\Bigl{(}\Phi^{-1}(u)-\Phi^{-1}(u_{i,k});0,b\Bigr{)}\cdot\phi\Bigl{(}\Phi^{-1}(v)-\Phi^{-1}(v_{j,k});0,b\Bigr{)}}{\phi\Bigl{(}\Phi^{-1}(u);0,b\Bigr{)}\cdot\phi\Bigl{(}\Phi^{-1}(v);0,b\Bigr{)}}\,$
(11)

where $\phi(.;0,b)$ is the pdf of a normal with mean $0$ and variance $b$,
$\Phi^{-1}$ is the inverse standard normal cdf. Samples
$\\{(u_{i,k},v_{j,k})\\}_{k=1}^{K}$ are easily obtained by iteratively fitting
KDE pair copulas on observed samples
$\\{(u_{1,k},\ldots,u_{d,k})\\}_{k=1}^{K}$ (Czado, 2019).

#### Choice of Hyperparameters.

We begin by choosing a Cauchy distribution as the initial predictive
distribution $P_{i}^{(0)}$ and the corresponding density $p_{i}^{(0)}\ \forall
i=1,\ldots,d$, with details provided in Appendix E. The hyperparameter $\rho$
for the R-BP recursion for each marginal predictive is chosen in an automatic
data-dependent manner, by minimising the energy score (Pacchiardi et al.,
2024) between observations and samples from our model conditional on previous
observed data. As the R-BP is sensitive to the ordering of the data, we follow
Fong et al. (2023); Ghalebikesabi et al. (2022) by averaging the resulting
R-BP marginal over 10 permutations of the data (see Tokdar et al. (2009);
Dixit and Martin (2019) for a discussion regarding the need of averaging over
permutations). We assume a the same variance parameter $b$ for all the KDE
pair copula estimators in the simplified vine and is chosen using 10-fold
cross-validation, in a data-dependent manner minimizing the same energy score.
The assumption of a common bandwidth $b$ is motivated by mapping all pair
copulas to a Gaussian space, which results in a common distance used on the
latent pair copula spaces. Another hyperparameter is the specific vine
decomposition (the grouping of dimensions in (10), see Appendix A.3) for which
we use a regular vine structure (Panagiotelis et al., 2017), selecting the
best pair-copula decomposition with a modified Bayesian Information Criterion
suited for vines (Nagler et al., 2019).

#### Benefits of the energy score

To note, unlike previous R-BP works (Hahn et al., 2018; Fong et al., 2023;
Ghalebikesabi et al., 2022) here we minimize the energy score to choose the
hyper-parameters rather than the log-score, both of which are strictly proper
scoring rules (see Appendix C and (Pacchiardi et al., 2024)) and define
statistical divergences. This choice was motivated by the robustness
properties of the energy score (Pacchiardi et al., 2021).

Further this choice also provides us some significant computational gain. We
only need samples from $P_{i}^{(n)}$ to get an unbiased estimate of the energy
score, which can be done very efficiently through inverse probability
sampling, hence halving the time to tune hyper-parameters compared to using
the log-score, which also requires densities $p_{i}^{(n)}$. Due to the
computational independence along each dimensions ensuing from Sklar’s copula
decomposition, this selection of $\rho$ (and model fitting) can further be
parallelised across dimensions (and permutations of the data) to exploit the
power of distributed computing. Similarly, the cross-validation of the copula
estimation can also be performed in parallel.

#### QB-Vine for Regression/classification tasks

Our framework can accommodate regression and classification tasks in addition
to density estimation, by rewriting the conditional density as following:

$p(y|\mathbf{x})=\frac{p(y,\mathbf{x})}{p(\mathbf{x})}=\frac{p_{y}(y)\cdot\prod_{i=1}^{d}\Bigl{\\{}p_{i}(x_{i})\Bigr{\\}}\cdot\mathbf{c}(y,x_{1},\ldots,x_{d})}{\prod_{i=1}^{d}\Bigl{\\{}p_{i}(x_{i})\Bigr{\\}}\cdot\mathbf{c}(x_{1},\ldots,x_{d})}=\frac{\mathbf{c}(y,x_{1},\ldots,x_{d})\cdot
p_{y}(y)}{\mathbf{c}(x_{1},\ldots,x_{d})}.$ (12)

The estimation of (12) is comprised of estimating the $d+1$ marginals for
$y,x_{1},\ldots,x_{d}$ and the two copulas $\mathbf{c}(y,\mathbf{x})$ and
$\mathbf{c}(\mathbf{x})$. We specify separate $\rho$ across marginal densities
as well as separate bandwidths for copulas to be estimated with the methods
described above. We note that for a copula decomposition to be unique, we
require that the marginals involved be continuous. This assumption is violated
in the case of classification for binary outcomes $y$. As such, we make use of
an approximation that transforms $y$ to a continuous scale by setting negative
examples to $-10$, and positive examples to $10$ and adds standard Gaussian
noise to all examples, breaking ties on a distribution scale (a common
approach taken in similar contexts (Journel and Xu, 1994; Liu, 2019; Harris et
al., 2022)). The rationale behind this approximation is that by setting the
two classes far apart on a marginal scale, we ensure no overlaps occur,
thereby maintaining a clear cut-off between the classes on the distribution
scale. Indeed, the separating boundary between the classes on a distribution
scale will be the percentile $q=\frac{T_{0}}{T_{1}+T_{0}}$ where $T_{0}$ and
$T_{1}$ are the numbers of negative and positive samples, respectively, in the
training set. Consequently, we create different clusters in the copula
$[0,1]^{d}$ hypercube according to the separation on the distributional scale,
facilitating the identification of patterns in the data. We note that other
approaches exist in the literature (Chen, 2016; Chen and Hanson, 2017; Czado
and Nagler, 2022) for classification with copulas which our framework can be
extended to.

#### Approximation error of Quasi-Bayesian Vine

To quantify the approximation power of the QB-Vine, we provide the following
stochastic boundedness (Bishop et al., 2007) result for univariate R-BP
distributions with respect to the limiting _martingale posterior_
$P^{(\infty)}$. To the best of our knowledge, this is the first convergence
rate result for the R-BP.

######  Lemma 3.2.

(R-BP predictive distribution convergence) The error of the distribution
function $P^{(n)}(x)$ in (4) is stochastically bounded with

$\sup_{x\in\mathcal{X}}\left|P^{(\infty)}(x)-P^{(n)}(x)\right|=\mathcal{O}_{p}\left(n^{-1/2}\right).$

###### Proof.

We provide a proof in Appendix D.1. ∎

In comparison to univariate Kernel Density Estimation with a mean-square
optimal bandwidth $b_{n}=\mathcal{O}(n^{-1/5})$, which convergence at a rate
$\mathcal{O}_{a.s.}(n^{-2/5}\sqrt{\text{ln}(n)})$, we can see that the
marginal R-BP has a better rate with sample size. In what follows, we assume
that the true copula is a simplified vine of which we know the decomposition
(a standard assumption in the vine copula litterature (Nagler and Czado, 2016;
Czado and Nagler, 2022; Spanhel and Kurz, 2019)). We strengthen marginal
guarantees with the theory on vine copulas to obtain the following convergence
result for the estimate of the copula density. In the statement of the
theorem, we consider marginal distributions $\\{P_{i}^{(\infty)}\\}_{i=1}^{d}$
and $\\{P_{i}^{(n)}\\}_{i=1}^{d}$ are implicitly applied to $\mathbf{x}$ for
respective copulas.

######  Theorem 3.3.

(Convergence of Quasi-Bayesian Vine) Assuming a correctly identified
simplified vine structure for $\mathbf{c}^{(\infty)}(\mathbf{u})$, and using
univariate R-BP marginal distributions with a simplified vine copula, the
copula estimator error is stochastically bounded
$\forall\,\mathbf{x}\in\mathbb{R}^{d}$ with

$|\mathbf{c}^{(\infty)}(\mathbf{x})-\mathbf{c}^{(n)}(\mathbf{x})|=\mathcal{O}_{p}(n^{-r})$
(13)

where $n^{-r}$ is the convergence rate of the KDE pair-copula.

###### Proof.

We provide a proof in Appendix D.2. ∎

For a bivariate KDE pair-copula estimator with optimal bandwidth
$b_{n}=\mathcal{O}\left(n^{-1/6}\right)$, we obtain $n^{-r}=n^{-1/3}$ (Nagler
and Czado, 2016). From Stone (1980), we note the optimal convergence rate of a
nonparametric estimator is $n^{-p/(2p+d)}$ where $p$ is the number of times
the estimator is differentiable. Therefore, as $d$ increases, we expect large
benefits from using a vine copula decomposition for the QB-Vine. When the
simplifying assumption does not hold, the vine copula converges to a partial
vine approximation of the true copula, as defined in Spanhel and Kurz (2019).
Together, these two results guarantee accurate samples from
$\mathbf{P}^{(n)}(\mathbf{x})$ by inverse probability sampling arguments. By
Theorem 3.3, the copula $\mathbf{c}^{(n)}$ ensures samples
$\mathbf{u}=(u_{1},\ldots,u_{d})$ on the $[0,1]^{d}$ hypercube have a
dependence structure representative of the data. Then, marginal distributions
$\\{P_{i}^{(n)}\\}_{i=1}^{d}$ recover dependent samples
$\mathbf{x}\in\mathbb{R}^{d}$ by evaluating the inverse of the distribution at
$u_{i}$ dimension-wise.

##  Related Work

Our method shares similarities with existing work on QB predictive density
estimation with analytical forms. The pivotal works of Newton et al. (1998);
Newton (2002) and the ensuing Predictive Recursion (PR) (Ghosh and Tokdar,
2006; Martin and Ghosh, 2008; Martin and Tokdar, 2009; Tokdar et al., 2009;
Martin, 2012; Ghosal and Van der Vaart, 2017; Martin, 2021) propose a
recursive solution to the same problem but are restrained to low dimensional
settings due to the numerical integration of a normalising constant over a
space scaling with $d$. A sequential importance sampling strategy for PR is
proposed in Dixit and Martin (2023) termed as PRticle Filter. The R-BP of Hahn
et al. (2018) and the multivariate extensions in Fong et al. (2023);
Ghalebikesabi et al. (2022) also have a recursive form driven by bivariate
copula updates. In the multivariate case, imposing assumptions on the kernel
structure leads to a conditional factorisation of the joint predictive which
recovers bivariate copula updates. In Ghalebikesabi et al. (2022), an
autoregressive Bayesian predictive (AR-BP) is used, where the dependence is
captured by dimension-wise similarity functions modelled with kernels or deep
autoregressive networks. The former relies on assumptions that might be too
simplistic to capture complex data while the latter loses the appeal of a
data-efficient predictive like in the R-BP. The Quasi-Bayesian Vine retains
the advantages of the bivariate copula-based recursion for marginal
predictives and circumvents the need for assumptions on the DPMM kernel. We
achieve this via approximating the high-dimensional dependency through a
simplified vine copula which is highly flexible and does not use a deep
network to preserve data-efficiency, all the while maintaining an analytical
expression. A relevant benchmark are the NFs of Papamakarios et al. (2017);
Durkan et al. (2019) with analytical densities with a solid performance across
data types and tasks.

##  Experiments

Table 1: Average log predictive score (lower is better) with error bars corresponding to two standard deviations over five runs for density estimation on datasets analysed by Ghalebikesabi et al. (2022). We note that as dimension increases, the QB-Vine outperforms all benchmarks. | WINE | BREAST | PARKIN | IONO | BOSTON
---|---|---|---|---|---
$n$/$d$ | 89/12 | 97/14 | 97/16 | 175/30 | 506/13
KDE | ${13.69_{\pm 0.00}}$ | ${10.45_{\pm 0.24}}$ | ${12.83_{\pm 0.27}}$ | ${32.06_{\pm 0.00}}$ | ${8.34_{\pm 0.00}}$
DPMM (Diag) | ${17.46_{\pm 0.60}}$ | ${16.26_{\pm 0.71}}$ | ${22.28_{\pm 0.66}}$ | ${35.30_{\pm 1.28}}$ | ${7.64_{\pm 0.09}}$
DPMM (Full) | ${32.88_{\pm 0.82}}$ | ${26.67_{\pm 1.32}}$ | ${39.95_{\pm 1.56}}$ | ${86.18_{\pm 10.22}}$ | ${9.45_{\pm 0.43}}$
MAF | ${39.60_{\pm 1.41}}$ | ${10.13_{\pm 0.40}}$ | ${11.76_{\pm 0.45}}$ | ${140.09_{\pm 4.03}}$ | ${56.01_{\pm 27.74}}$
RQ-NSF | ${38.34_{\pm 0.63}}$ | ${26.41_{\pm 0.57}}$ | ${31.26_{\pm 0.31}}$ | ${54.49_{\pm 0.65}}$ | ${-2.20_{\pm 0.11}}$
PRticle Filter | ${23.89_{\pm 0.93}}$ | ${25.98_{\pm 1.06}}$ | ${34.79_{\pm 3.95}}$ | ${79.22_{\pm 9.87}}$ | ${27.18_{\pm 3.12}}$
R-BP | ${13.57_{\pm 0.04}}$ | ${7.45_{\pm 0.02}}$ | ${9.15_{\pm 0.04}}$ | ${21.15_{\pm 0.04}}$ | ${4.56_{\pm 0.04}}$
Rd-BP | ${13.32_{\pm 0.01}}$ | ${6.12_{\pm 0.05}}$ | ${7.52_{\pm 0.05}}$ | ${19.82_{\pm 0.08}}$ | ${-13.50_{\pm 0.59}}$
AR-BP | ${13.45_{\pm 0.05}}$ | ${6.18_{\pm 0.05}}$ | ${8.29_{\pm 0.11}}$ | ${17.16_{\pm 0.25}}$ | ${-0.45_{\pm 0.77}}$
ARd-BP | $\mathbf{13.22_{\pm 0.04}}$ | ${6.11_{\pm 0.04}}$ | ${7.21_{\pm 0.12}}$ | ${16.48_{\pm 0.26}}$ | ${-14.75_{\pm 0.89}}$
ARnet-BP | ${14.41_{\pm 0.11}}$ | ${6.87_{\pm 0.23}}$ | ${8.29_{\pm 0.17}}$ | ${15.32_{\pm 0.35}}$ | ${-5.71_{\pm 0.62}}$
QB-Vine | ${13.76_{\pm 0.13}}$ | $\mathbf{4.67_{\pm 0.31}}$ | $\mathbf{4.93_{\pm 0.20}}$ | $\mathbf{-16.08_{\pm 2.12}}$ | $\mathbf{-31.04_{\pm 1.02}}$

In this section, we compare our QB-Vine model against competing methods
supporting density evaluation with a closed-form expression. Further details
on the experiments are included in Appendix E, with an ablation study and
comparison with sample size and dimension. Code is included in the
supplementary material.

#### Density estimation

We evaluate the QB-Vine on density estimation benchmark UCI datasets (Asuncion
and Newman, 2007) with small sample sizes ranging from $89$ to $506$ and
dimensionality varying from $12$ to $30$, adding results for the QB-Vine and
PRticle Filter to the experiments of Ghalebikesabi et al. (2022). We report
the log predictive score
LPS$=\frac{1}{n_{test}}\sum_{k}^{n_{test}}-\mathbf{p}^{(n_{train})}(\mathbf{x}_{k})$
on a held-out test dataset of size $n_{test}$ comprised of half the samples
with the other half used for training, averaging results over five runs with
random partitions each time. We compare the QB-Vine against the following
models: Kernel Density Estimation (Parzen, 1962), DPMM (Rasmussen, 1999) with
a diagonal (Diag) and full (Full) covariance matrix for each mixture
component, MAF (Papamakarios et al., 2017), RQ-NSF (Durkan et al., 2019) as
well as the closely related PRticle Filter (Dixit and Martin, 2023), R-BP
(Fong et al., 2023) and AR-BP (Ghalebikesabi et al., 2022). For the last two
Bayesian predictive models, we add a subscript $d$ to indicate that the $\rho$
hyperparameter possibly differs across dimensions, and the $net$ suffix
indicates a network-based selection of $\rho$ for dimensions. We observe in
Table 1 that our QB-Vine method comfortably outperforms all competitors as the
dimension increases, matching the performance of other Bayesian predictive
models for the lower dimensional WINE dataset. Our method’s relative
performance increases with the dimension of the data, particularly achieving a
much smaller LPS for IONO - the dataset with the largest dimensions and a
relatively small sample size. We accredit this performance to the copula
decomposition as that is our main distinguishing factor from the Bayesian
predictive models.

#### High-dimensional image dataset

Figure 2: Density estimation on the Digits data ($n=1797,d=64$) with reduced
training sizes for the QB-Vine against other models fitted on the full
training set. The QB-Vine achieves competitive performance for training sizes
as little as $n=50$ and outperforms all competitors once $n>200$.

We further evaluate the QB-Vine on the digits dataset ($n=$1797, $d$=64) as a
high-dimensional example with a relatively low sample size. The high contrast
between $n$ and $d$ makes the problem suited for assessing the data efficiency
and convergence of the QB-Vine. We compare with the two NF models as their
high model capacity is a good fit for image data, as well as all the Bayesian
predictive methods of Section 5, from the study of Ghalebikesabi et al.
(2022). We report the average LPS in bits per dimension (bpd) with standard
errors over five runs with random partitions, using half the sample size to
train models and the other half to evaluate the LPS. Additionally, we report
the average LPS of the QB-Vine, obtained in the same way except for the
training set size being reduced (to $30,50,100,200,300,400,500$). Figure 2
depicts the QB-Vine’s performance for different-sized training sets. When
trained on the full train set, the QB-Vine outperforms all competitors by a
considerable margin. Furthermore, our method is competitive with as little as
50 training samples and outperforms all benchmarks past a training size of
200, demonstrating its data-efficiency and convergence speed. A complete
numerical table is reported in Appendix E.

Table 2: Average LPS (lower is better) with error bars corresponding to two standard deviations over five runs for supervised tasks analysed by Ghalebikesabi et al. (2022). The QB-Vine performs favourably against benchmarks, with relative performance improving as samples per dimension decrease. | Regression |  | Classification
---|---|---|---
| BOSTON | CONCR | DIAB |  | IONO | PARKIN
$n$/$d$ | 506/13 | 1,030/8 | 442/10 |  | 351/33 | 195/22
Linear | ${0.87_{\pm 0.03}}$ | ${0.99_{\pm 0.01}}$ | ${1.07_{\pm 0.01}}$ |  | ${0.33_{\pm 0.01}}$ | ${0.38_{\pm 0.01}}$
GP | ${0.42_{\pm 0.08}}$ | ${0.36_{\pm 0.02}}$ | ${1.06_{\pm 0.02}}$ |  | ${0.30_{\pm 0.02}}$ | ${0.42_{\pm 0.02}}$
MLP | ${1.42_{\pm 1.01}}$ | ${2.01_{\pm 0.98}}$ | ${3.32_{\pm 4.05}}$ |  | ${0.26_{\pm 0.05}}$ | ${0.31_{\pm 0.02}}$
R-BP | ${0.76_{\pm 0.09}}$ | ${0.87_{\pm 0.03}}$ | ${1.05_{\pm 0.03}}$ |  | ${0.26_{\pm 0.01}}$ | ${0.37_{\pm 0.01}}$
Rd-BP | ${0.40_{\pm 0.03}}$ | ${0.42_{\pm 0.00}}$ | ${1.00_{\pm 0.02}}$ |  | ${0.34_{\pm 0.02}}$ | ${0.27_{\pm 0.03}}$
AR-BP | ${0.52_{\pm 0.13}}$ | ${0.42_{\pm 0.01}}$ | ${1.06_{\pm 0.02}}$ |  | ${0.21_{\pm 0.02}}$ | ${0.29_{\pm 0.02}}$
ARd-BP | ${0.37_{\pm 0.10}}$ | ${0.39_{\pm 0.01}}$ | ${0.99_{\pm 0.02}}$ |  | ${0.20_{\pm 0.02}}$ | ${0.28_{\pm 0.03}}$
ARnet-BP | ${0.45_{\pm 0.11}}$ | $\mathbf{-0.03_{\pm 0.00}}$ | ${1.41_{\pm 0.07}}$ |  | ${0.24_{\pm 0.04}}$ | ${0.26_{\pm 0.04}}$
QB-Vine | $\mathbf{-0.81_{\pm 1.26}}$ | ${0.54_{\pm 0.34}}$ | $\mathbf{0.87_{\pm 0.20}}$ |  | $\mathbf{-1.85_{\pm 1.16}}$ | $\mathbf{-0.76_{\pm 0.28}}$

#### Regression and classification

We further demonstrate our method’s effectiveness on supervised learning
tasks, with three datasets for regression and two datasets for classification,
adding to the study of Ghalebikesabi et al. (2022). For classification, we
transform the binary values to continuous ones to preserve copula assumptions,
as detailed in Section 3. We report the conditional LPS
$=\frac{1}{n_{test}}\sum_{k}^{n_{test}}-\mathbf{p}^{(n_{train})}(y_{k}|\mathbf{x}_{k})$

over a test set of size $n_{test}$ made up of half the samples with the
conditional estimator trained on the other half of the data. We compare our
model against a Gaussian Process (Williams and Rasmussen, 1995), a linear
Bayesian model (Linear) (Minka, 2000), a one-hidden-layer multilayer
perceptron (MLP), as well the R-BP and AR-BP variants for supervised tasks
(Fong et al., 2023; Ghalebikesabi et al., 2022). The QB-Vine outperforms
competing methods on all datasets except CONCR. We believe the lower
performance on CONCR is due to the high number of samples relative to the
dimension, preventing our approach from fully exploiting the vine copula
decomposition. Once again, the performance of the QB-Vine more clearly exceeds
that of competitors as dimensions increase. The QB-Vine has higher standard
errors than other methods (except MLP), which we posit is the consequence of
our conditional estimator in Section 3 being defined as a ratio, inflating the
variation in the LPS. However, we highlight that an overly precise inference
is more misleading/dangerous than an overly uncertain one.

##  Discussion

We introduced the Quasi-Bayesian Vine, a joint Bayesian predictive density
estimator with an analytical form and easy to sample. This extends the
existing works on Quasi-Bayesian predictive densities, by using Sklar’s
theorem to decompose the predictive density into predictive marginals and a
copula to model the high-dimensional dependency. This decomposition also
provides us a way for two-part estimation procedure, employing recursive
density estimation as in the Quasi-Bayesian works for the marginals and
fitting a simplified vine copula for the dependence, resulting in a
convergence rate independent of dimension for certain joint densities. We
empirically demonstrate the advantage of QB-Vine on a range of datasets
compared to other benchmark methods, showing excellent modeling capabilities
in large dimensions with only a few training data.

However, there is potential for further improvements. The non-uniqueness of
the vine decomposition results in a search over a vast model space during
estimation. In addition, the hyperparameter selection of the KDE pair copulas
could be ameliorated. Finally, the main assumption is the use of a simplified
vine copula which is an approximation to the true distribution. From a
practical point of view, however, it offers fast computations and outperforms
competitors as shown in experiments. Other, future directions of this work
include incorporation of more effective copula models, or copulas
accommodating different dependence structures (Vatter and Nagler, 2018; Nagler
et al., 2022; Huk et al., 2023a). Another exciting direction are developments
of new recursive Quasi-Bayes methods that can be merged into a Quasi-Bayesian
Vine model (Garelli et al., 2024).

## References

  * Alquier et al. [2022] Pierre Alquier, Badr-Eddine Chérief-Abdellatif, Alexis Derumigny, and Jean-David Fermanian.  Estimation of copulas via maximum mean discrepancy.  _Journal of the American Statistical Association_ , pages 1–16, 2022.
  * Arbel et al. [2021] Michael Arbel, Alex Matthews, and Arnaud Doucet.  Annealed flow transport Monte Carlo.  In _International Conference on Machine Learning_ , pages 318–330. PMLR, 2021.
  * Ashok et al. [2024] Arjun Ashok, Étienne Marcotte, Valentina Zantedeschi, Nicolas Chapados, and Alexandre Drouin.  TACTis-2: Better, faster, simpler attentional copulas for multivariate time series.  In _The Twelfth International Conference on Learning Representations_ , 2024.  URL https://openreview.net/forum?id=xtOydkE1Ku.
  * Asuncion and Newman [2007] Arthur Asuncion and David Newman.  UCI machine learning repository, 2007.
  * Bedford and Cooke [2001] Tim Bedford and Roger M Cooke.  Probability density decomposition for conditionally dependent random variables modeled by vines.  _Annals of Mathematics and Artificial intelligence_ , 32:245–268, 2001.
  * Bellemare et al. [2017] Marc G Bellemare, Ivo Danihelka, Will Dabney, Shakir Mohamed, Balaji Lakshminarayanan, Stephan Hoyer, and Rémi Munos.  The Cramer distance as a solution to biased Wasserstein gradients.  _arXiv preprint arXiv:1705.10743_ , 2017.
  * Berti et al. [2004] Patrizia Berti, Luca Pratelli, and Pietro Rigo.  Limit theorems for a class of identically distributed random variables.  2004\.
  * Berti et al. [2006] Patrizia Berti, Luca Pratelli, and Pietro Rigo.  Almost sure weak convergence of random probability measures.  _Stochastics and Stochastics Reports_ , 78(2):91–97, 2006.
  * Berti et al. [2023] Patrizia Berti, Emanuela Dreassi, Fabrizio Leisen, Luca Pratelli, and Pietro Rigo.  A probabilistic view on predictive constructions for Bayesian learning.  _Statistical Science_ , 1(1):1–15, 2023.
  * Beskos et al. [2013] Alexandros Beskos, Natesh Pillai, Gareth Roberts, Jesus-Maria Sanz-Serna, and Andrew Stuart.  Optimal tuning of the hybrid Monte Carlo algorithm.  2013\.
  * Bishop et al. [2007] Yvonne M Bishop, Stephen E Fienberg, and Paul W Holland.  _Discrete multivariate analysis: Theory and practice_.  Springer Science & Business Media, 2007.
  * Charpentier et al. [2007] Arthur Charpentier, Jean-David Fermanian, and Olivier Scaillet.  The estimation of copulas: Theory and practice.  _Copulas: From theory to application in finance_ , 35, 2007.
  * Chen [2016] Yuhui Chen.  A copula-based supervised learning classification for continuous and discrete data.  _Journal of Data Science_ , 14(4):769–782, 2016.
  * Chen and Hanson [2017] Yuhui Chen and Timothy Hanson.  Copula regression models for discrete and mixed bivariate responses.  _Journal of Statistical Theory and Practice_ , 11:515–530, 2017.
  * Clayton [1978] David G Clayton.  A model for association in bivariate life tables and its application in epidemiological studies of familial tendency in chronic disease incidence.  _Biometrika_ , 65(1):141–151, 1978.
  * Czado [2019] Claudia Czado.  Analyzing dependent data with vine copulas.  _Lecture Notes in Statistics, Springer_ , 222, 2019.
  * Czado and Nagler [2022] Claudia Czado and Thomas Nagler.  Vine copula based modeling.  _Annual Review of Statistics and Its Application_ , 9:453–477, 2022.
  * Dawid et al. [2016] A Philip Dawid, Monica Musio, and Laura Ventura.  Minimum scoring rule inference.  _Scandinavian Journal of Statistics_ , 43(1):123–138, 2016.
  * De Finetti [1937] Bruno De Finetti.  La prévision: ses lois logiques, ses sources subjectives.  In _Annales de l’institut Henri Poincaré_ , volume 7, pages 1–68, 1937.
  * Dixit and Martin [2019] Vaidehi Dixit and Ryan Martin.  Permutation-based uncertainty quantification about a mixing distribution.  _arXiv preprint arXiv:1906.05349_ , 2019.
  * Dixit and Martin [2023] Vaidehi Dixit and Ryan Martin.  A PRticle filter algorithm for nonparametric estimation of multivariate mixing distributions.  _Statistics and Computing_ , 33(4):1–14, 2023.
  * Durkan et al. [2019] Conor Durkan, Artur Bekasov, Iain Murray, and George Papamakarios.  Neural spline flows.  _Advances in neural information processing systems_ , 32, 2019.
  * Elidan [2013] Gal Elidan.  Copulas in machine learning.  In _Copulae in Mathematical and Quantitative Finance: Proceedings of the Workshop Held in Cracow, 10-11 July 2012_ , pages 39–60. Springer, 2013.
  * Fischer et al. [2009] Matthias Fischer, Christian Köck, Stephan Schlüter, and Florian Weigert.  An empirical analysis of multivariate copula models.  _Quantitative Finance_ , 9(7):839–854, 2009.
  * Fong et al. [2023] Edwin Fong, Chris Holmes, and Stephen G Walker.  Martingale posterior distributions.  _Journal of the Royal Statistical Society Series B: Statistical Methodology_ , 85(5):1357–1391, 2023.
  * Fortini and Petrone [2020] Sandra Fortini and Sonia Petrone.  Quasi-Bayes properties of a procedure for sequential learning in mixture models.  _Journal of the Royal Statistical Society Series B: Statistical Methodology_ , 82(4):1087–1114, 2020.
  * Fortini and Petrone [2023] Sandra Fortini and Sonia Petrone.  Prediction-based uncertainty quantification for exchangeable sequences.  _Philosophical Transactions of the Royal Society A_ , 381(2247):20220142, 2023.
  * Frazier et al. [2024] David T. Frazier, Jeremias Knoblauch, and Christopher Drovandi.  The impact of loss estimation on Gibbs measures, 2024.
  * Garelli et al. [2024] Samuele Garelli, Fabrizio Leisen, Luca Pratelli, and Pietro Rigo.  Asymptotics of predictive distributions driven by sample means and variances, 2024.
  * Geenens et al. [2017] Gery Geenens, Arthur Charpentier, and Davy Paindaveine.  Probit transformation for nonparametric kernel estimation of the copula density.  2017\.
  * Ghalebikesabi et al. [2022] Sahra Ghalebikesabi, Chris Holmes, Edwin Fong, and Brieuc Lehmann.  Density estimation with autoregressive Bayesian predictives.  _arXiv preprint arXiv:2206.06462_ , 2022.
  * Ghosal and Van der Vaart [2017] Subhashis Ghosal and Aad Van der Vaart.  _Fundamentals of nonparametric Bayesian inference_ , volume 44.  Cambridge University Press, 2017.
  * Ghosh and Tokdar [2006] Jayanta K Ghosh and Surya T Tokdar.  Convergence and consistency of Newton’s algorithm for estimating mixing distribution.  In _Frontiers in statistics_ , pages 429–443. World Scientific, 2006.
  * Gneiting and Raftery [2007] Tilmann Gneiting and Adrian E Raftery.  Strictly proper scoring rules, prediction, and estimation.  _Journal of the American statistical Association_ , 102(477):359–378, 2007.
  * Gretton et al. [2006] Arthur Gretton, Karsten Borgwardt, Malte Rasch, Bernhard Schölkopf, and Alex Smola.  A kernel method for the two-sample-problem.  _Advances in neural information processing systems_ , 19, 2006.
  * Größer and Okhrin [2022] Joshua Größer and Ostap Okhrin.  Copulae: An overview and recent developments.  _Wiley Interdisciplinary Reviews: Computational Statistics_ , 14(3):e1557, 2022.
  * Gu et al. [2015] Shixiang Shane Gu, Zoubin Ghahramani, and Richard E Turner.  Neural adaptive sequential Monte Carlo.  _Advances in neural information processing systems_ , 28, 2015.
  * Haff et al. [2010] Ingrid Hobæk Haff, Kjersti Aas, and Arnoldo Frigessi.  On the simplified pair-copula construction—simply useful or too simplistic?  _Journal of Multivariate Analysis_ , 101(5):1296–1310, 2010.
  * Hahn et al. [2018] P Richard Hahn, Ryan Martin, and Stephen G Walker.  On recursive Bayesian predictive distributions.  _Journal of the American Statistical Association_ , 113(523):1085–1093, 2018.
  * Harris et al. [2022] Lucy Harris, Andrew TT McRae, Matthew Chantry, Peter D Dueben, and Tim N Palmer.  A generative deep learning approach to stochastic downscaling of precipitation forecasts.  _Journal of Advances in Modeling Earth Systems_ , 14(10):e2022MS003120, 2022.
  * Hewitt and Savage [1955] Edwin Hewitt and Leonard J Savage.  Symmetric measures on Cartesian products.  _Transactions of the American Mathematical Society_ , 80(2):470–501, 1955.
  * Hjort et al. [2010] Nils Lid Hjort, Chris Holmes, Peter Müller, and Stephen G Walker.  _Bayesian nonparametrics_ , volume 28.  Cambridge University Press, 2010.
  * Hofert et al. [2013] Marius Hofert, Martin Mächler, and Alexander J McNeil.  Archimedean copulas in high dimensions: Estimators and numerical challenges motivated by financial applications.  _Journal de la Société Française de Statistique_ , 154(1):25–63, 2013.
  * Hofert et al. [2018] Marius Hofert, Ivan Kojadinovic, Martin Mächler, and Jun Yan.  _Elements of copula modeling with R_.  Springer, 2018.
  * Hoffman et al. [2019] Matthew Hoffman, Pavel Sountsov, Joshua V Dillon, Ian Langmore, Dustin Tran, and Srinivas Vasudevan.  Neutralizing bad geometry in Hamiltonian Monte Carlo using neural transport.  _arXiv preprint arXiv:1903.03704_ , 2019.
  * Holmes and Walker [2023] Chris Holmes and Stephen G Walker.  Statistical inference with exchangeability and martingales.  _Philosophical Transactions of the Royal Society A_ , 381(2247):20220143, 2023.
  * Huk et al. [2023a] David Huk, Rilwan A Adewoyin, and Ritabrata Dutta.  Probabilistic rainfall downscaling: Joint generalized neural models with censored spatial Gaussian copula.  _arXiv preprint arXiv:2308.09827_ , 2023a.
  * Huk et al. [2023b] David Huk, Lorenzo Pacchiardi, Ritabrata Dutta, and Mark Steel.  David Huk, Lorenzo Pacchiardi, Ritabrata Dutta and Mark Steel’s contribution to the discussion of ‘Martingale posterior distributions’ by Fong, Holmes and Walker.  _Journal of the Royal Statistical Society Series B: Statistical Methodology_ , 85(5):1405–1406, 2023b.
  * Janke et al. [2021] Tim Janke, Mohamed Ghanmi, and Florian Steinke.  Implicit generative copulas.  _Advances in Neural Information Processing Systems_ , 34:26028–26039, 2021.
  * Joe [2014] Harry Joe.  _Dependence modeling with copulas_.  CRC press, 2014.
  * Journel and Xu [1994] André G Journel and Wenlong Xu.  Posterior identification of histograms conditional to local data.  _Mathematical Geology_ , 26:323–359, 1994.
  * Kamthe et al. [2021] Sanket Kamthe, Samuel Assefa, and Marc Deisenroth.  Copula flows for synthetic data generation.  _arXiv preprint arXiv:2101.00598_ , 2021.
  * Kantorovich [1960] Leonid V Kantorovich.  Mathematical methods of organizing and planning production.  _Management science_ , 6(4):366–422, 1960.
  * Kingma and Welling [2013] Diederik P Kingma and Max Welling.  Auto-encoding variational Bayes.  _arXiv preprint arXiv:1312.6114_ , 2013.
  * Ling et al. [2020] Chun Kai Ling, Fei Fang, and J Zico Kolter.  Deep Archimedean copulas.  _Advances in Neural Information Processing Systems_ , 33:1535–1545, 2020.
  * Liu [2019] Weiwei Liu.  Copula multi-label learning.  _Advances in Neural Information Processing Systems_ , 32, 2019.
  * Lueckmann et al. [2021] Jan-Matthis Lueckmann, Jan Boelts, David Greenberg, Pedro Goncalves, and Jakob Macke.  Benchmarking simulation-based inference.  In _International conference on artificial intelligence and statistics_ , pages 343–351. PMLR, 2021.
  * Martin [2012] Ryan Martin.  Convergence rate for predictive recursion estimation of finite mixtures.  _Statistics & Probability Letters_, 82(2):378–384, 2012.
  * Martin [2021] Ryan Martin.  A survey of nonparametric mixing density estimation via the predictive recursion algorithm.  _Sankhya B_ , 83:97–121, 2021.
  * Martin and Ghosh [2008] Ryan Martin and Jayanta K Ghosh.  Stochastic approximation and Newton’s estimate of a mixing distribution.  _Statistical Science_ , pages 365–382, 2008.
  * Martin and Tokdar [2009] Ryan Martin and Surya T Tokdar.  Asymptotic properties of predictive recursion: robustness and rate of convergence.  2009\.
  * Matthews et al. [2022] Alex Matthews, Michael Arbel, Danilo Jimenez Rezende, and Arnaud Doucet.  Continual repeated annealed flow transport Monte Carlo.  In _International Conference on Machine Learning_ , pages 15196–15219. PMLR, 2022.
  * Minka [2000] Thomas Minka.  Bayesian linear regression.  Technical report, Citeseer, 2000.
  * Nagler and Czado [2016] Thomas Nagler and Claudia Czado.  Evading the curse of dimensionality in nonparametric density estimation with simplified vine copulas.  _Journal of Multivariate Analysis_ , 151:69–89, 2016.
  * Nagler et al. [2017] Thomas Nagler, Christian Schellhase, and Claudia Czado.  Nonparametric estimation of simplified vine copula models: comparison of methods.  _Dependence Modeling_ , 5(1):99–120, 2017.
  * Nagler et al. [2019] Thomas Nagler, Christian Bumann, and Claudia Czado.  Model selection in sparse high-dimensional vine copula models with an application to portfolio risk.  _Journal of Multivariate Analysis_ , 172:180–192, 2019.
  * Nagler et al. [2022] Thomas Nagler, Daniel Krüger, and Aleksey Min.  Stationary vine copula models for multivariate time series.  _Journal of Econometrics_ , 227(2):305–324, 2022.
  * Newton [2002] Michael A Newton.  On a nonparametric recursive estimator of the mixing distribution.  _Sankhyā: The Indian Journal of Statistics, Series A_ , pages 306–322, 2002.
  * Newton et al. [1998] Michael A Newton, Fernando A Quintana, and Yunlei Zhang.  Nonparametric Bayes methods using predictive updating.  In _Practical nonparametric and semiparametric Bayesian statistics_ , pages 45–61. Springer, 1998.
  * Pacchiardi and Dutta [2022] Lorenzo Pacchiardi and Ritabrata Dutta.  Score matched neural exponential families for likelihood-free inference.  _The Journal of Machine Learning Research_ , 23(1):1745–1815, 2022.
  * Pacchiardi et al. [2021] Lorenzo Pacchiardi, Sherman Khoo, and Ritabrata Dutta.  Generalized Bayesian likelihood-free inference.  _arXiv preprint arXiv:2104.03889_ , 2021.
  * Pacchiardi et al. [2024] Lorenzo Pacchiardi, Rilwan A Adewoyin, Peter Dueben, and Ritabrata Dutta.  Probabilistic forecasting with generative networks via scoring rule minimization.  _Journal of Machine Learning Research_ , 25(45):1–64, 2024.
  * Panagiotelis et al. [2017] Anastasios Panagiotelis, Claudia Czado, Harry Joe, and Jakob Stöber.  Model selection for discrete regular vine copulas.  _Computational Statistics & Data Analysis_, 106:138–152, 2017.
  * Papamakarios et al. [2017] George Papamakarios, Theo Pavlakou, and Iain Murray.  Masked autoregressive flow for density estimation.  _Advances in neural information processing systems_ , 30, 2017.
  * Papamakarios et al. [2019] George Papamakarios, David Sterratt, and Iain Murray.  Sequential neural likelihood: Fast likelihood-free inference with autoregressive flows.  In _The 22nd international conference on artificial intelligence and statistics_ , pages 837–848. PMLR, 2019.
  * Papamakarios et al. [2021] George Papamakarios, Eric Nalisnick, Danilo Jimenez Rezende, Shakir Mohamed, and Balaji Lakshminarayanan.  Normalizing flows for probabilistic modeling and inference.  _Journal of Machine Learning Research_ , 22(57):1–64, 2021.
  * Parzen [1962] Emanuel Parzen.  On estimation of a probability density function and mode.  _The annals of mathematical statistics_ , 33(3):1065–1076, 1962.
  * Rasmussen [1999] Carl Rasmussen.  The infinite Gaussian mixture model.  _Advances in Neural Information Processing Systems_ , 12, 1999.
  * Rezende et al. [2014] Danilo Jimenez Rezende, Shakir Mohamed, and Daan Wierstra.  Stochastic backpropagation and approximate inference in deep generative models.  In _International conference on machine learning_ , pages 1278–1286. PMLR, 2014.
  * Roberts and Rosenthal [1998] Gareth O Roberts and Jeffrey S Rosenthal.  Optimal scaling of discrete approximations to Langevin diffusions.  _Journal of the Royal Statistical Society: Series B (Statistical Methodology)_ , 60(1):255–268, 1998.
  * Roberts et al. [1997] GO Roberts, A Gelman, and WR Gilks.  Weak convergence and optimal scaling of random walk Metropolis algorithms.  _The annals of applied probability_ , 7(1):110–120, 1997.
  * Schepsmeier et al. [2015] Ulf Schepsmeier, Jakob Stoeber, Eike Christian Brechmann, Benedikt Graeler, Thomas Nagler, Tobias Erhardt, Carlos Almeida, Aleksey Min, Claudia Czado, Mathias Hofmann, et al.  Package ‘vinecopula’.  _R package version_ , 2(5), 2015.
  * Sethuraman [1994] Jayaram Sethuraman.  A constructive definition of Dirichlet priors.  _Statistica sinica_ , pages 639–650, 1994.
  * Si et al. [2023] Phillip Si, Zeyi Chen, Subham Sekhar Sahoo, Yair Schiff, and Volodymyr Kuleshov.  Semi-autoregressive energy flows: exploring likelihood-free training of normalizing flows.  In _International Conference on Machine Learning_ , pages 31732–31753. PMLR, 2023.
  * Sklar [1959] M Sklar.  Fonctions de répartition à n dimensions et leurs marges.  In _Annales de l’ISUP_ , volume 8, pages 229–231, 1959.
  * Song et al. [2017] Jiaming Song, Shengjia Zhao, and Stefano Ermon.  A-nice-mc: Adversarial training for mcmc.  _Advances in neural information processing systems_ , 30, 2017.
  * Spanhel and Kurz [2019] Fabian Spanhel and Malte S Kurz.  Simplified vine copula models: Approximations based on the simplifying assumption.  _Electronic Journal of Statistics_ , 13:1254–1291, 2019.
  * Sriperumbudur et al. [2009] Bharath K Sriperumbudur, Arthur Gretton, Kenji Fukumizu, Gert RG Lanckriet, and Bernhard Schölkopf.  A note on integral probability metrics and $\varphi$-divergences.  _arXiv preprint arXiv:0901.2698_ , 2009.
  * Stone [1980] Charles J Stone.  Optimal rates of convergence for nonparametric estimators.  _The annals of Statistics_ , pages 1348–1360, 1980.
  * Székely and Rizzo [2005] Gábor J Székely and Maria L Rizzo.  A new test for multivariate normality.  _Journal of Multivariate Analysis_ , 93(1):58–80, 2005.
  * Tokdar et al. [2009] Surya T Tokdar, Ryan Martin, and Jayanta K Ghosh.  Consistency of a recursive estimate of mixing distributions.  _The Annals of Statistics_ , pages 2502–2522, 2009.
  * Vatter and Nagler [2018] Thibault Vatter and Thomas Nagler.  Generalized additive models for pair-copula constructions.  _Journal of Computational and Graphical Statistics_ , 27(4):715–727, 2018.
  * Williams and Rasmussen [1995] Christopher Williams and Carl Rasmussen.  Gaussian processes for regression.  _Advances in neural information processing systems_ , 8, 1995.

##  Appendix A Copulas

Copulas are a widely adopted tool in statistics and machine learning for
modelling densities, permitting the construction of joint densities in a two-
step estimation process. One firstly estimates the marginal densities as if
they were independent of each other, and secondly, models the _copula_ which
accounts for the dependence between dimensions. This approach is motivated
through Sklar’s theorem:

######  Theorem A.1 (Sklar Sklar [1959]).

Let $\mathbf{P}$ be an $d$-dimensional distribution function with continuous
marginal distributions $P_{1},P_{2},\ldots,P_{d}$. Then there exists a copula
distribution $\mathbf{C}$ such that for all
$\mathbf{x}=(x_{1},x_{2},\ldots,x_{d})\in\mathbb{R}^{d}$:

$\mathbf{P}(x_{1},\ldots,x_{d})=\mathbf{C}(P_{1}(x_{1}),\ldots,P_{n}(x_{d}))$
(14)

And if a probability density function (pdf) is available:

$\mathbf{p}(x_{1},\ldots,x_{d})=p_{1}(x_{1})\cdot\ldots\cdot
p_{d}(x_{d})\cdot\mathbf{c}(P_{1}(x_{1}),\ldots,P_{d}(x_{d}))$ (15)

where $p_{1}(x_{1}),\ldots,p_{d}(x_{d})$ are the marginal pdfs, and
$\mathbf{c}(P_{1}(x_{1}),\ldots,P_{d}(x_{d}))$ is the copula pdf.

If the marginal distributions are absolutely continuous, the copula is unique.
Consequently, one can decompose the estimation problem of learning
$p(\mathbf{x})$ by first learning all the marginals $\\{p_{i}\\}_{i=1}^{d}$,
and in a second step learning an appropriate copula model
$c(u_{1},\ldots,u_{d})$, where $u_{i}:=P_{i}(x_{i}),i\in\\{1,\ldots d\\}$ are
the images of the $x_{i}$ under the cdf of each dimension. By applying cdf
transformations marginally, the copula is agnostic to the differences between
dimensions such as axis scaling, and purely focuses upon capturing the
dependence structure among them.

Most parametric copula models are only suited for two-dimensional dependence
modelling and greatly suffer from the curse of dimensionality [Fischer et al.,
2009, Hofert et al., 2013]. The Gaussian copula is a popular parametric choice
as it is well-studied and can be fitted quickly even to moderate dimensions
[Hofert et al., 2018]. However, it lacks the desired flexibility to capture
more complex dependencies involving multiple dimensions. Among nonparametric
copulas, Kernel Density Estimator (KDE) copulas [Geenens et al., 2017] are
commonly used. They apply an inverse Gaussian distribution to the observed
$\\{u_{i}\\}_{i=1}^{d}$ to map them to a latent space and perform regular KDE
on the latent density. However, this KDE copula method suffers from the poor
scaling of KDE estimators in higher dimensions. Finally, deep learning copula
models remain a nascent line of research and typically are more
computationally expensive and sample-dependent due to their reliance on large
models, with only a handful of candidate solutions such as Ling et al. [2020],
Kamthe et al. [2021], Janke et al. [2021], Ashok et al. [2024]. As such,
current copula models are mostly limited to low to medium-dimensional
modelling [Joe, 2014].

###  A.1 Gaussian copula

A popular parametric copula model is the Gaussian copula. It assumes that the
dependence between dimensions is identical to that of a Gaussian distribution
with mean $\mathbf{0}$ and covariance matrix $\Sigma$:

$c(u_{1},\ldots,u_{d})=\frac{\mathcal{N}_{d}(\Phi^{-1}(u_{1}),\ldots,\Phi^{-1}(u_{d});\mathbf{0},\Sigma)}{\prod_{i=1}^{d}\mathcal{N}(\Phi^{-1}(u_{i});0,1)}.$
(16)

As such its only parameters are the off-diagonal entries of the covariance
matrix $\Sigma$. In the case of $d=2$, there is a single parameter to estimate
for the Gaussian copula.

###  A.2 Gaussian Kernel Density Estimator copulas

As the marginal distribution of a copula is uniform in $[0,1]$, the support of
the copula estimator is restricted to $[0,1]^{d}$ and must satisfy the uniform
marginal condition. It is fairly difficult to build such an estimator that
fulfills both desiderata with high expressivity. Gaussian Kernel Density
Estimation for copulas [Charpentier et al., 2007, Geenens et al., 2017] is a
popular approach for such nonparametric copula models, which models the copula
on a latent Gaussian space. We explain the approach in the following.

By the inverse sampling theorem, for any $\mathbb{R}$-valued continuous random
variable $X$, applying the corresponding cumulative distribution function $F$
to $X$ results in $F(X)$ being uniformly distributed in $[0,1]$. Thus, we can
transform a uniform random variable into any continuous distribution using its
inverse cumulative distribution function. In the copula estimation stage,
samples from the copula already have uniform marginal distributions, meaning
we can apply any inverse distribution $F^{-1}$ to each marginal sample value
and obtain a corresponding latent marginal distribution. If $F^{-1}$ is the
inverse standard normal distribution, then the latent distribution for each
marginal will be normal and $\mathbb{R}$-valued with no uniformity
restrictions.

Gaussian KDE for copulas applies inverse standard normal distributions to each
marginal, resulting in a latent representation of the samples on a Gaussian
space. As such, one can employ regular Gaussian KDE to estimate the copula
density on this latent space. In the case of two-dimensional copulas, the
ensuing estimator has the following expression:

$\hat{c}(u,v)\,=\,\frac{\sum_{k=1}^{K}\phi\Bigl{(}\Phi^{-1}(u)-\Phi^{-1}(u_{k});0,b\Bigr{)}\cdot\phi\Bigl{(}\Phi^{-1}(v)-\Phi^{-1}(v_{k});0,b\Bigr{)}}{\phi\Bigl{(}\Phi^{-1}(u);0,b\Bigr{)}\cdot\phi\Bigl{(}\Phi^{-1}(v);0,b\Bigr{)}}\,,$
(17)

where $(u,v)$ and $(u_{k},v_{k})$ are both in
$[0,1]^{2}$,$\\{(u_{k},v_{k})\\}_{k=1}^{K}$ are observed copula samples, and
$\Phi$ and $\phi$ are respectively the Gaussian distribution and density with
mean $0$ and variance $b>0$.

###  A.3 Vine copulas

A vine copula is an efficient dependency modelling method which decomposes
$d$-dimensional copula estimation into $d(d-1)/2$ bivariate copula estimation
via structured conditioning [Bedford and Cooke, 2001]. Here we illustrate the
decomposition by a $3$-dimensional copula density $c_{U,V,W}$ for the random
vector ($U,\,V,\,W$):

$c_{U,V,W}(u,v,w)\,=\,c_{U,V}(u,v)\cdot c_{V,W}(v,w)\cdot
c_{U,W|V}\bigg{(}C_{U|V}(u|\,v),\,C_{W|V}(w|\,v)\bigg{|}\,v\bigg{)}\,,$ (18)

where

  * •

$C_{U,V}$ is the copula of $(U,\,V)$, $C_{V,W}$ is the copula of $(V,\,W)$;

  * •

$C_{U,W|V=v}$ is the copula of $(C_{U|V}(U|\,v),C_{W|V}(W|\,v))$ conditional
on $V=v$;

  * •

$C_{U|V}$ is the conditional distribution of $C_{U,V}$ on $V$, and $C_{W|V}$
is the conditional distribution of $C_{V,W}$ on $V$.

Generally, the distribution of $C_{U,W|V=v}$ will change with different values
$v$ and this will make the model relatively complex. Therefore, it is common
to use the simplifying assumption by ignoring the conditioning of pair copulas
and simply model them as unconditional bivariate densities:

$c_{U,V,W}(u,v,w)\,=\,c_{U,V}(u,v)\cdot c_{V,W}(v,w)\cdot
c_{U,W|V}\bigg{(}C_{U|V}(u|\,v),\,C_{W|V}(w|\,v)\bigg{)}\,.$ (19)

The rationale of the simplified assumption is studied in Haff et al. [2010].
In this paper, we mainly focus on the regular vine copula (R-vine) with a
simplified assumption. The construction of an R-vine copula has two basic
ingredients: (1) a valid tree structure of all random variables, (2) the
choice of family for bivariate copulas.

Before we introduce the tree structure which is valid to construct a R-vine
copula, let us first rigorously define the random variable of a copula we want
to estimate. Suppose $\mathbf{U=(U_{1},U_{2},...,U_{d})}$ is a $d$-dimensional
random variable which is distributed as a copula $C$, then we have that each
marginal random variable $U_{i}$ is uniformly distributed in $[0,1]$-scale.
For notational simplicity, we will use the index of a random variable as its
notation instead. Denote
$\mathcal{T}=\\{\mathcal{T}_{1},...,\mathcal{T}_{d-1}\\}$ as a sequence of
$(d-1)$ trees in terms of $\mathcal{T}_{i}=(\mathcal{V}_{i},\mathcal{E}_{i})$.
Here we use a set of two nodes to represent the corresponding edge in
$\mathcal{T}_{i}$, i.e., $e=(a,b)$ if node $a$ and $b$ in $\mathcal{E}_{i}$
are linked. To construct a valid R-vine copula, $\mathcal{T}$ satisfies the
following conditions:

  * •

$\mathcal{T}_{1}$ is a tree with a set of edges $\mathcal{E}_{1}$ and a set of
nodes $\mathcal{V}_{1}=\\{1,...,d\\}$;

  * •

$\mathcal{T}_{i}$ is a tree with a set of edges $\mathcal{E}_{i}$ and a set of
nodes $\mathcal{V}_{i}=\mathcal{V}(\mathcal{E}_{i-1})$ for $i=2,3,...,d-1$,
where $\mathcal{V}(E)$ denoted that the pair of nodes which are linked by an
edge in $E$ is treated as a new node.

  * •

If two nodes in $\mathcal{E}_{i+1}$ are linked by an edge in
$\mathcal{V}_{i+1}$, then they must be linked to one common node in
$\mathcal{T}_{i}$.

For $\forall\,e=(a,b)\in\mathcal{T}_{i}$ with $i\geq 2$, we define

$c(e)=a\cap b,\quad\quad a(e)=a\setminus c(e),\quad\quad b(e)=b\setminus
c(e)\,.$ (20)

Finally, we rigorously define the R-vine copula as follows.

######  Definition A.2 (Regular Vine Copulas).

A $d$-dimensional copula $C$ is a regular vine copula if there exists a tuple
$(\mathcal{T},\mathcal{C})$ such that

  * •

$\mathcal{T}$ is a regular vine tree sequence with $(d-1)$ trees;

  * •

$\mathcal{C}=\\{C_{e}:e\in\mathcal{E}_{i},i\in[d-1],C_{e}\text{ is a bivariate
copula}\\}$ is a family of bivariate copulas for each edge;

  * •

For $\forall e\in\mathcal{E}_{i}$ and $\forall i\in[d-1]$, then $C_{e}$ is
corresponding to the copula of $\bigg{(}a(e),b(e)\bigg{)}\bigg{|}\,c(e)$.

Therefore, the density function of R-vine $C$ can be expressed as

$\begin{split}&c_{(\mathcal{T},\mathcal{C})}(u_{1},...,u_{d})\\\
&=\prod_{i=1}^{d-1}\prod_{e\in\mathcal{E}_{i}}c_{(a(e),b(e))|c(e)}(C_{a(e)|c(e)}(v_{a(e)}|v_{c(e)}),C_{b(e)|c(e)}(v_{b(e)}|v_{c(e)}))\,.\end{split}$
(21)

Here we illustrate an R-vine copula density in five dimensions where we use
different colors for the pair copulas corresponding to each of the $d-1=4$
trees.

$\begin{split}c_{\mathcal{T},\mathcal{C}}(u_{1},u_{2},u_{3},u_{4},u_{5})=&c(u_{1},u_{2})\cdot
c(u_{1},u_{3})\cdot c(u_{2},u_{4})\cdot c(u_{3},u_{5})\\\
&{\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}\cdot
c_{1,5|3}(u_{1|3},u_{5|3})\cdot c_{2,3|1}(u_{2|1},u_{3|1})\cdot
c_{1,4|2}(u_{1|2},u_{4|2})}\\\
&{\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}\cdot
c_{2,5|1,3}(u_{2|1,3},u_{5|1,3})\cdot c_{3,4|1,2}(u_{3|1,2},u_{4|1,2})}\\\
&{\color[rgb]{0,.5,.5}\definecolor[named]{pgfstrokecolor}{rgb}{0,.5,.5}\cdot
c_{4,5|1,2,3}(u_{4|1,2,3},u_{5|1,2,3})}\,.\end{split}$ (22)

Specifying the tree structure for an R-vine decomposition is essential and
plays an important role in pair-copula estimation through conditioning. An
excellent overview is given in Czado [2019]. Notably, an R-vine decomposition
is not unique for a given joint copula pdf. An appealing tree selection
algorithm is proposed in Nagler et al. [2019], where authors derive a modified
BIC criterion that prioritizes sparse trees while being consistent when the
dimension $d$ grows at most at the rate of $\sqrt{n}$ where $n$ is the sample
size.

##  Appendix B Martingale Posterior Distributions

Here we explain _martingale posterior_ distributions as a justification of the
Bayesian approach through a focus on prediction.

###  B.1 The Bayesian Choice as a Consequence of Predictive Uncertainty

A common goal in statistics is the inference of a parameter or quantity
$\theta$ by analysing data in the form of observations
$(x_{1},\ldots,x_{n}),\,n\in\mathbb{N}$. The rationale for learning from data
is that each observation provides information about the underlying process and
parameter $\theta$, without which statistical analysis would be redundant.
Indeed, consider a decision-maker basing their decision on their belief about
$\theta$ (in an i.i.d. setting). Having observed data $(x_{1},\ldots,x_{n})$
and given the opportunity to observe an additional data point $x_{n+1}$, they
would be assumed to accept, as they could refine their beliefs on $\theta$.
This process of updating one’s beliefs based on data is at the core of the
Bayesian approach. Equipped with an initial guess about the parameter of
interest $\theta$ captured by the prior $\pi(\theta)$, the goal is the
inference about the distribution of the parameter given observed data
$(x_{1},\ldots,x_{n})$ which is encoded in the posterior density
$\pi^{(n)}\left(\theta|(x_{1},\ldots,x_{n})\right)$.

For the decision-maker to refuse, it would mean that the additional
observation has no significant effect on their belief. This implies, as
identified by Fong et al. [2023] and Holmes and Walker [2023], that there is a
point where additional observations $(x_{n+1},\ldots,x_{N})$ provide no
benefit to the knowledge update of $\theta$. Inspecting the Bayesian posterior
in terms of observed data $x_{1:n}=(x_{1},\ldots,x_{n})$ and possible future
observations to be made $x_{n+1:\infty}$, written as

$\pi^{(n)}\left(\theta|x_{1:n}\right)=\int\pi^{(n)}\left(\theta,x_{n+1:\infty}|x_{1:n}\right)~{}dx_{n+1:\infty},$
(23)

one can expand the right-hand integrand by including the predictive density
$p$ for future data, obtaining

$\pi^{(n)}\left(\theta|x_{1:n}\right)=\int\underbrace{\pi_{\infty}\left(\theta|x_{1:\infty}\right)}_{\text{Bayes
estimate}}\cdot\underbrace{p(x_{n+1:\infty}|x_{1:n})}_{\text{predictive
density}}~{}dx_{n+1:\infty}.$ (24)

Having rewritten the posterior density in this way, it becomes apparent that
the uncertainty in the value of $\theta$, which is given by the Bayes
estimate, is a consequence of the uncertainty surrounding the imputation of
missing observations $x_{n+1:\infty}$ through the predictive. With this
insight, unlike the traditional Bayesian construction of prior-likelihood,
Fong et al. [2023] proceed to replace the predictive density $p$ with a
predictive mechanism using all available data to impute the missing
$x_{n+1:\infty}$ (or at least impute $x_{n+1:N}$ for a sufficiently large
$N$), and replace the Bayes estimate $\pi^{(\infty)}$ with an appropriate
functional of the complete data $x_{1:\infty}$. This predictive mechanism used
to impute the unobserved data $x_{n+1:\infty}$ given observed $x_{1:n}$
directly leads to the _martingale posterior_ as the limiting distribution of
functionals when the unobserved data has been imputed.

###  B.2 Constructing a martingale posterior for the DPMM

As explained in Section 2, given a sequence of observations $(x_{i})_{i\geq
1}$ from the data generating process, the Bayesian approach wants to obtain an
infinite sequence of copula densities $(c_{i})_{i\geq 1}$ to recursively
update the Bayesian predictive distribution. However, finding such an infinite
sequence of copulas isn’t feasible in practice and the ensuing copula family
is determined by the prior-likelihood pair one implicitly chooses. For
example, from Hahn et al. [2018], we have the following copula sequences:

  * •

If we select the likelihood and prior as $l(y;\theta)=\theta e^{-\theta y}$
and $\pi(\theta)=e^{-\theta}$, then

$c_{n}(u,v_{n})=\frac{(n+1)[(1-u)^{-1-1/n}(1-v_{n})^{-1-1/n}]}{n[(1-u)^{-1/n}+(1-v_{n})^{-1/n}-1]^{n+2}}\,,$
(25)

which is a sequence of Clayton copula [Clayton, 1978] with parameter $n^{-1}$.

  * •

If we select the likelihood and prior as $l(y|\theta)=\phi(y;\theta,1)$ and
$\pi(\theta)=\phi(\theta;0,\tau^{-1})$, then

$c_{n}(u,v_{n})=\frac{\phi_{2}(\Phi^{-1}(u),\Phi^{-1}(v_{n});\bm{0},\bm{\Sigma}_{\rho_{n}})}{\phi(\Phi^{-1}(u);0,1)\phi(\Phi^{-1}(v_{n});0,1)}\,,$
(26)

where $\phi_{2}$ is the pdf of a bivariate normal distribution, $\phi$ is the
pdf of normal distribution, and

$\bm{\Sigma}_{\rho_{n}}=\begin{bmatrix}1&\rho_{n}\\\
\rho_{n}&1\end{bmatrix}\,,\quad\quad\rho_{n}=(n+\tau)^{-1}\,.$ (27)

Therefore, we obtain a sequence of Gaussian copulas with parameters
$\\{\rho_{n}\\}_{n\geq 1}$. Further, if we additionally put a conjugate prior
on the variance parameter $\sigma^{2}$ of the likelihood
$l(y|\theta)=\phi(y;\theta,\sigma^{2})$, then we will recover a sequence of
Student-t copulas.

The issue of model mis-specification naturally arises here due to the
selection of prior-likelihood pair. Fortunately, here we can employ the DPMM
which has relatively high flexibility due to its nonparametric nature.For the
DPMM, suppose the first observation $x_{1}$ arrives, then we can obtain the
first predictive via the updating kernel

$\displaystyle k_{1}(x,x_{1})$
$\displaystyle=\frac{\mathbb{E}\left[f_{G}(x)\,f_{G}(x_{1})\right]}{\mathbb{E}[f_{G}(x)]\cdot\mathbb{E}[f_{G}(x_{1})]}\,,$
(28)

where $p_{0}(x)=\mathbb{E}[f_{G}(x)]$. Then, we can derive the numerator as

$\displaystyle\mathbb{E}\left[f_{G}(x)\,f_{G}(x_{1})\right]$ (29)
$\displaystyle=$
$\displaystyle\mathbb{E}\left[\sum_{j=1}^{\infty}\sum_{k=1}^{\infty}w_{j}\,w_{k}\,\phi\left(x;\theta_{j},1\right)\phi\left(x_{1};\theta_{k},1\right)\right]$
$\displaystyle=$
$\displaystyle\left(1-\mathbb{E}\left[\sum_{i=1}^{\infty}w_{i}^{2}\right]\right)\mathbb{E}_{G_{0}}\left[\phi(x;\theta,1)\right]\,\mathbb{E}_{G_{0}}\left[\phi(x_{1};\theta,1)\right]$
$\displaystyle+$
$\displaystyle\mathbb{E}\left[\sum_{i=1}^{\infty}w_{i}^{2}\right]\,\mathbb{E}_{G_{0}}\left[\phi(x;\theta,1)\,\phi(x_{1};\theta,1)\right]$
$\displaystyle=$
$\displaystyle\left(1-\mathbb{E}\left[\sum_{i=1}^{\infty}w_{i}^{2}\right]\right)p_{0}(x)\,p_{0}(x_{1})+\mathbb{E}\left[\sum_{i=1}^{\infty}w_{i}^{2}\right]\,\mathbb{E}_{G_{0}}\left[\phi(x;\theta,1)\,\phi(x_{1};\theta,1)\right]\,.$

The first equality follows from the stick-breaking representation [Sethuraman,
1994] of the $\operatorname{DP}$, we can formulate $G$ as

$G(\,\cdot\,)\,=\,\sum_{i=1}^{\infty}w_{i}\,\delta_{\theta_{i}}(\,\cdot\,)$
(30)

where

$w_{i}=v_{i}\prod_{j<i}(1-v_{j}),\quad\quad v_{i}\stackrel{{\scriptstyle
i.i.d.}}{{\sim}}\operatorname{Beta}(1,a),\quad\quad\theta_{i}\stackrel{{\scriptstyle
i.i.d.}}{{\sim}}G_{0}.$ (31)

Then, the second equality follows from the condition that
$\sum_{i=1}^{\infty}w_{i}=1$ almost surely. Denote
$\alpha=\mathbb{E}\left[\sum_{i=1}^{\infty}w_{i}^{2}\right]$, then we can
write the updating kernel as

$\displaystyle k_{1}(x,x_{1})$
$\displaystyle=(1-\alpha)+\alpha\cdot\frac{\mathbb{E}_{G_{0}}\left[\phi(x;\theta,1)\,\phi(x_{1};\theta,1)\right]}{p_{0}(x)\,p_{0}(x_{1})}$
(32)
$\displaystyle=(1-\alpha)+\alpha\cdot\frac{\phi_{2}(\Phi^{-1}(u),\Phi^{-1}(v_{1});\bm{0},\Sigma_{\rho})}{\phi(\Phi^{-1}(u);0,1)\,\phi(\Phi^{-1}(v_{1});0,1)}$
$\displaystyle=(1-\alpha)+\alpha\cdot
c_{\rho}(\Phi^{-1}(u),\Phi^{-1}(v_{1}))\,,$

where

  * •

covariance matrix $\Sigma_{\rho}=\begin{bmatrix}1&\rho\\\
\rho&1\end{bmatrix}$,

  * •

$\Phi^{-1}$ is the inverse cdf of standard normal distribution,

  * •

$u=\mathbb{P}_{0}(x)$, $v_{1}=\mathbb{P}_{0}(x_{1})$,

  * •

$p_{0}(\,\cdot\,)=\phi(\,\cdot\,;0,1+\tau^{-1})$,

  * •

$c_{\rho}$ is a bivariate Gaussian copula density function with the
correlation parameter $\rho=1/(1+\tau)$.

Therefore, we can see that the copula $c_{1}(u,v_{1})$ for the first updating
step is a mixture of independent copula and Gaussian copula. However, we will
lose the tractable form of the copula from the second updating step. Following
Hahn et al. [2018], instead of deriving the explicit rule for the sequence of
copula, we fix the correlation parameter $\rho$ and set $\alpha$ to be a
$(0,1)$-valued decreasing sequence $(\alpha_{i})_{i\geq 1}$.

##  Appendix C Strictly Proper Scoring Rules

In probabilistic machine learning, a Scoring Rule (SR) measures the
appropriateness of a distribution $\mathbb{P}$ in modelling an observation
$\mathbf{x}\in\mathcal{X}$ through a score written as
$S(\mathbb{P},\mathbf{x})$. For (hyper) parameter estimation, suppose we aim
to model the underlying distribution of an observation $\mathbf{x}$ using a
family of distributions $\mathbb{P}_{\mathbf{\theta}}$ parametrized by
$\mathbf{\theta}$, then we can use a SR $S$ to select the suitable
$\mathbf{\theta}$. If we assume that $\mathbf{x}\sim\mathbb{Q}$, then we can
obtain the expected SR $\mathcal{S}$ via taking expectation w.r.t.
$\mathbf{x}$ as

$\mathcal{S}(\mathbb{P},\mathbb{Q})\,=\,\mathbb{E}_{\mathbf{x}\sim\mathbb{Q}}\Bigl{[}S(\mathbb{P},\mathbf{x})\Bigr{]}\,.$
(33)

Following Gneiting and Raftery [2007], we call $\mathcal{S}$ strictly proper
if $\mathcal{S}(\mathbb{P},\mathbb{Q})$ is minimized if and only if
$\mathbb{P}=\mathbb{Q}$, i.e. for $\forall\,\mathbb{P}\in\mathcal{P}$ with
$\mathbb{P}\neq\mathbb{Q}$ such that

$\mathcal{S}(\mathbb{Q},\mathbb{Q})<\mathcal{S}(\mathbb{P},\mathbb{Q})\,.$
(34)

Considering our settings of marginal distributions, we have observed data
$\mathbf{x}_{1:n}\stackrel{{\scriptstyle i.i.d.}}{{\sim}}\mathbb{P}^{*}$, and
we aim to model $\mathbb{P}^{*}$ using $\mathbb{P}_{\mathbf{\theta}}$. More
explicitly, we don’t need to use $\mathbb{P}_{\mathbf{\theta}}$ directly in
the expected SR, instead we normally use its probability density function or
samples from this distribution to evaluate. In general, we optimise

$\displaystyle\mathbf{\theta}^{*}$
$\displaystyle=\underset{\mathbf{\theta}\in\Theta}{\text{argmin}}\mathbb{E}_{\mathbf{\theta}\in\Theta}\mathcal{S}(\mathbb{P}_{\mathbf{\theta}},\mathbb{Q})$
$\displaystyle=\underset{\mathbf{\theta}\in\Theta}{\text{argmin}}\,\mathcal{S}(\mathbb{P}_{\mathbf{\theta}},\mathbb{Q}).$

Since we do not have the complete population of $\mathbb{Q}$, we use the
empirical SR $\hat{\mathcal{S}}$ instead, i.e.

$\hat{\mathbf{\theta}}^{*}\,=\,\underset{\mathbf{\theta}\in\Theta}{\text{argmin}}\mathbb{E}_{\mathbf{\theta}\in\Theta}\hat{\mathcal{S}}(\mathbb{P}_{\mathbf{\theta}},\mathbb{Q})\,,$
(35)

where
$\hat{\mathcal{S}}(\mathbb{P}_{\mathbf{\theta}},\mathbb{Q})=\frac{1}{n}\sum_{k=1}^{n}S(\mathbb{P}_{\mathbf{\theta}},\mathbf{x}_{k})$.
Under mild conditions, it can be proven that
$\hat{\mathbf{\theta}}^{*}\rightarrow\mathbf{\theta}^{*}$ asymptotically in
Dawid et al. [2016]. For any positive definite kernel
$k(\,\cdot\,,\,\cdot\,)$, the Kernel Score [Gneiting and Raftery, 2007] is
given by

$S_{k}(\mathbb{P}_{\mathbf{\theta},\mathbf{x}})\,=\,\mathbb{E}[k(\mathbf{Y},\mathbf{Y}^{\prime})]-2\cdot\mathbb{E}[k(\mathbf{Y},\mathbf{x})]\,,$
(36)

where $\mathbf{Y},\mathbf{Y}^{\prime}\sim\mathbb{P}_{\mathbf{\theta}}$. If we
set $k(\mathbf{x},\mathbf{y})=-||\mathbf{x}-\mathbf{y}||^{\beta}_{2}$, then we
obtain the energy score which is strictly proper if
$\mathbb{E}_{\mathbf{Y}\sim\mathbb{P}_{\mathbf{\theta}}}||\mathbf{Y}||^{\beta}_{2}<\infty$.
The energy Score is defined as

$S_{E}^{\beta}(\mathbb{P}_{\mathbf{\theta}},\mathbf{x})\,=\,2\cdot\mathbb{E}||\mathbf{Y}-\mathbf{x}||^{\beta}_{2}-\mathbb{E}||\mathbf{Y}-\mathbf{Y}^{\prime}||^{\beta}_{2}\,,\quad\text{for
}\beta\in(0,2]\,.$ (37)

Here $\mathbf{Y}and\mathbf{Y}^{\prime}$ are i.i.d. samples from
$\mathbb{P}_{\mathbf{\theta}}$. Practically, given finite samples
$\mathbf{y}_{1}(\mathbf{\theta}),...,\mathbf{y}_{m}(\mathbf{\theta})\stackrel{{\scriptstyle
i.i.d.}}{{\sim}}\mathbb{P}_{\mathbf{\theta}}$ where
$\mathbf{y}(\mathbf{\theta})$ denotes a sample from
$\mathbb{P}_{\mathbf{\theta}}$ which is a differentiable function of
$\mathbf{\theta}$, the unbiased estimate of the energy Score via Monte Carlo
approximation is

$\hat{S}_{E}^{\beta}(\mathbf{y}_{1:m}(\mathbf{\theta}),\mathbf{x})=\frac{2}{m}\sum_{j=1}^{m}||\mathbf{y}_{j}(\mathbf{\theta})-\mathbf{x}||^{\beta}_{2}-\frac{1}{m(m-1)}\sum_{k\neq
j}||\mathbf{y}_{j}(\mathbf{\theta})-\mathbf{y}_{k}(\mathbf{\theta})||^{\beta}_{2}\,.$
(38)

Similarly, we can derive the unbiased gradient of
$S_{E}^{\beta}(\mathbf{y}_{1:m}(\mathbf{\theta}),\mathbf{x})$ w.r.t.
$\mathbf{\theta}$ which is crucial for any gradient descent algorithm. For
notational convenience, we define

$g(\mathbf{Y},\mathbf{Y}^{\prime},\mathbf{x})=2\cdot||\mathbf{Y}-\mathbf{x}||^{\beta}_{2}-||\mathbf{Y}-\mathbf{Y}^{\prime}||^{\beta}_{2}\,,$
(39)

then
$S_{E}^{\beta}(\mathbb{P}_{\mathbf{\theta}},\mathbf{x})=\mathbb{E}_{\mathbf{Y},\mathbf{Y}^{\prime}\sim\mathbb{P}_{\mathbf{\theta}}}[g(\mathbf{Y},\mathbf{Y}^{\prime},\mathbf{x})]$.
Next, we have that

$\begin{split}&\nabla_{\mathbf{\theta}}S_{E}^{\beta}(\mathbb{P}_{\mathbf{\theta}},\mathbf{x})\\\
=&\nabla_{\mathbf{\theta}}\mathbb{E}_{\mathbf{Y},\mathbf{Y}^{\prime}\sim\mathbb{P}_{\mathbf{\theta}}}[g(\mathbf{Y},\mathbf{Y}^{\prime},\mathbf{x})]\\\
=&\mathbb{E}_{\mathbf{Y},\mathbf{Y}^{\prime}\sim\mathbb{P}_{\mathbf{\theta}}}\Bigl{[}\nabla_{\mathbf{\theta}}g(\mathbf{Y},\mathbf{Y}^{\prime},\mathbf{x})\Bigr{]}\\\
\simeq&\frac{1}{m(m-1)}\sum_{j=1}^{m}\sum_{k=1}^{m}\nabla_{\mathbf{\theta}}g(\mathbf{y}_{j}(\mathbf{\theta}),\mathbf{y}_{k}(\mathbf{\theta}),\mathbf{x})\cdot\delta_{\\{j\neq
k\\}}\\\
=&\frac{2}{m}\sum_{j=1}^{m}\nabla_{\mathbf{\theta}}||\mathbf{y}_{j}(\mathbf{\theta})-\mathbf{x}||^{\beta}_{2}-\frac{1}{m(m-1)}\sum_{k\neq
j}\nabla_{\mathbf{\theta}}||\mathbf{y}_{j}(\mathbf{\theta})-\mathbf{y}_{k}(\mathbf{\theta})||^{\beta}_{2}\\\
=&\widehat{\nabla}_{\mathbf{\theta}}S_{E}^{\beta}(\mathbf{y}_{1:m}(\mathbf{\theta}),\mathbf{x})\,.\end{split}$
(40)

Furthermore, if $\mathbb{P}_{\mathbf{\theta}}$ and $\mathbb{Q}$ both are
distributions of $\mathbb{R}$-valued random variables, then the energy Score
will be reduced to

$S_{E}^{\beta}(\mathbb{P}_{\mathbf{\theta}},x)\,=\,2\cdot\mathbb{E}|Y-x|^{\beta}-\mathbb{E}|Y-Y^{\prime}|^{\beta}\,,\quad\text{for
}\beta\in(0,2]\,,$ (41)

where $Y,Y^{\prime}\sim\mathbb{P}_{\mathbf{\theta}}$. Notice that this will
become the Continuous Ranked Probability Score (CRPS) [Székely and Rizzo,
2005] when $\beta=1$.

The energy Score is a strictly proper scoring rule [Gneiting and Raftery,
2007, Dawid et al., 2016, Pacchiardi and Dutta, 2022] and is a special
instance of the maximum mean discrepancy [Gretton et al., 2006] as well as a
statistical divergence. It has been used as an effective objective for copula
estimation [Janke et al., 2021, Alquier et al., 2022, Huk et al., 2023a] and
even for R-BP marginal predictives as shown in Huk et al. [2023b]. It has also
enjoyed success as a objective for Normalising Flows [Si et al., 2023] and
generative models [Pacchiardi et al., 2024]. The energy Score is the only
objective among Wasserstein $p$-metric [Kantorovich, 1960] that supports
unbiased gradient evaluations [Bellemare et al., 2017, Pacchiardi and Dutta,
2022], has a known optimisation-free solution and features a faster
convergence rate than similar integral probability metrics [Sriperumbudur et
al., 2009, Frazier et al., 2024].

##  Appendix D Proofs

###  D.1 Lemma 3.2

###### Proof.

We follow the Definition 14.4-3 of stochastic boundedness from Bishop et al.
[2007]. Begin by choosing $\delta\in(0,1)$. From Proposition 1 in Fong et al.
[2023], we have for $\epsilon>0$ and $M>n$, over the supremum of
$x\in\mathcal{X}$:

$\begin{split}&\mathbb{P}\left(\left|P^{(M)}(x)-P^{(n)}(x)\right|\geq\epsilon\right)\leq
2\exp\left(-\frac{\epsilon^{2}}{\frac{2\epsilon\alpha_{n+1}}{3}+\frac{1}{2}\sum_{i=n+1}^{M}\alpha_{i}^{2}}\right)\\\
\Leftrightarrow\;\lim_{M\rightarrow\infty}&\mathbb{P}\left(\left|P^{(M)}(x)-P^{(n)}(x)\right|\leq\epsilon\right)\geq\lim_{M\rightarrow\infty}1-2\exp\left(-\frac{\epsilon^{2}}{\frac{2\epsilon\alpha_{n+1}}{3}+\frac{1}{2}\sum_{i=n+1}^{M}\alpha_{i}^{2}}\right).\end{split}$
(42)

Next, we choose a value $\epsilon_{n}$ (where the subscript shows that this
quantity is dependent on $n$) to have the appropriate probability on the
right-hand side by enforcing:

$\begin{split}&\delta=2\exp\left(-\frac{\epsilon_{n}^{2}}{\frac{2\epsilon_{n}\alpha_{n+1}}{3}+\frac{1}{2}\sum_{i=n+1}^{\infty}\alpha_{i}^{2}}\right)\\\
\Leftrightarrow\;&0=\epsilon_{n}^{2}+\log\left(\frac{\delta}{2}\right)\frac{2\alpha_{n+1}}{3}\epsilon_{n}+\frac{1}{2}\log\left(\frac{\delta}{2}\right)\sum_{i=n+1}^{\infty}\alpha_{i}^{2}\end{split}$
(43)

with solution

$\epsilon_{n}=\frac{-\log\left(\frac{\delta}{2}\right)\frac{2\alpha_{n+1}}{3}+\sqrt{\left[\log\left(\frac{\delta}{2}\right)\frac{2\alpha_{n+1}}{3}\right]^{2}-2\log\left(\frac{\delta}{2}\right)\sum_{i=n+1}^{\infty}\alpha_{i}^{2}}}{2}.$
(44)

Due to $\delta\in(0,1)$ and $\alpha_{i}>0\,\forall i$, we have:

$|\epsilon_{n}|=n^{-1/2}\frac{-\log\left(\frac{\delta}{2}\right)n^{1/2}\frac{2\alpha_{n+1}}{3}+\sqrt{\left[\log\left(\frac{\delta}{2}\right)n^{1/2}\frac{2\alpha_{n+1}}{3}\right]^{2}-2\log\left(\frac{\delta}{2}\right)n\sum_{i=n+1}^{\infty}\alpha_{i}^{2}}}{2}.$
(45)

With the choice of $\alpha_{i}=(2-\frac{1}{i})(\frac{1}{i+1})$, we see that
$\lim_{n\rightarrow\infty}n^{-a}\alpha_{i}$ is bounded for all powers
$a\geq-1$. As such, both
$-\log\left(\frac{\delta}{2}\right)n^{1/2}\frac{2\alpha_{n+1}}{3}$ and
$\left[\log\left(\frac{\delta}{2}\right)n^{1/2}\frac{2\alpha_{n+1}}{3}\right]^{2}$
will safely be bounded for large enough $n$. Similarly, due to the choice of
$\alpha_{i}$, we have
$\sum_{i=n+1}^{\infty}(\alpha_{i})^{2}=\mathcal{O}(n^{-1})$, meaning
$n^{-a}\sum_{i=n+1}^{\infty}(\alpha_{i})^{2}$ will be bounded for large enough
$n$ as long as $a\geq-1$. Hence
$-2\log\left(\frac{\delta}{2}\right)n\sum_{i=n+1}^{\infty}\alpha_{i}^{2}$ will
also be bounded for large enough $n$.

Consequently, for our choice of $\delta$, there exists a finite $K>0$ and a
finite $N>0$ such that:

$\sup_{x\in\mathcal{X}}\mathbb{P}\left(\left|\frac{P^{(\infty)}(x)-P^{(n)}(x)}{n^{-1/2}}\right|\leq
K\right)\geq\ 1-\delta\quad\forall n>N.$ (46)

∎

###  D.2 Theorem 3.3

We prove the statement for general densities, which then naturally extends to
predictive densities. The proof of this result is largely an adaptation of the
simplified vine copula convergence result (Theorem 1 in Nagler and Czado
[2016]) but where no rate on marginal densities is required. As such, our
proof shares an identical approach until the last part, where we deviate. We
have a weaker result for the convergence of distributions and yet show in what
follows that convergence of our copula estimator can be obtained even without
convergence guarantees on marginal densities.

#### Notation used in the proof

We follow vine copula notation from Appendix A.3 and use superscripts with
parenthesis to now differentiate between samples instead of predictive steps,
following notational conventions of the literature. We define $h$-functions as
the conditional distribution functions for pair copulas

$h_{j_{e}\mid\ell_{e};D_{e}^{\prime}}(u\mid
v):=\int_{0}^{u}c_{j_{e},\ell_{e};D_{e}^{\prime}}(s,v)ds,\quad\text{ for
}(u,v)\in[0,1]^{2}.$ (47)

Further, we refer to the true unobserved samples of pair-copulas as

$U_{j_{e}\mid D_{e}}^{(i)}:=F_{je}\mid
D_{e}\left(X_{j_{e}}^{(i)}\mid\bm{X}_{D_{e}}^{(i)}\right),\quad U_{k_{e}\mid
D_{e}}^{(i)}:=F_{k_{e}\mid
D_{e}}\left(X_{k_{e}}^{(i)}\mid\bm{X}_{D_{e}}^{(i)}\right),$ (48)

for $i=1,\ldots,n$. We also denote estimators and quantities obtained by
application to these unobserved samples with a bar superscript, for example:

$\bar{c}_{j_{e},k_{e};D_{e}}(u,v):=\bar{c}_{j_{e},k_{e};D_{e}}\left(u,v,U_{j_{e}\mid
D_{e}}^{(1)},\ldots,U_{k_{e}\mid D_{e}}^{(n)}\right).$ (49)

Finally, denote with a hat superscript all quantities and estimators obtained
by using $\hat{U_{l}}^{(i)}:=\hat{F_{l}}(X_{l}^{(i)})$ instead of the true
unobserved samples used in Equation (48).

#### Assumptions

For completeness, we state assumptions about the marginal distribution
estimator $\hat{P}$ as well as bivariate copula estimators $\hat{c}$, even
though our practical choices from the main text respect these. We begin by
stating the assumption about our marginal distribution estimator denoted
$\hat{P}$:

  * •

A1: The marginal distribution function estimator has the following convergence
rate:

$\sup_{x\in\mathcal{X}}\left|\hat{P}(x)-P(x)\right|=\mathcal{O}_{p}\left(n^{-1/2}\right).$
(50)

Next, we state our assumptions about the pair copula estimator for
completeness:

  * •

A2: For all $e\in E_{m},m=1,\ldots,d-1$, with $-r$ the convergence rate of a
bivariate KDE copula estimator, it holds:
(a) for all $(u,v)\in(0,1)^{2}$,

$\bar{c}_{j_{e},k_{e};D_{e}}(u,v)-c_{j_{e},k_{e};D_{e}}(u,v)=O_{p}\left(n^{-r}\right),$
(51)

(b) for every $\delta\in(0,0.5]$,

$\displaystyle\sup_{(u,v)\in[\delta,1-\delta]^{2}}\left|\bar{h}_{j_{e}\mid
k_{e};D_{e}}(u\mid v)-h_{j_{e}\mid k_{e};D_{e}}(u\mid v)\right|$
$\displaystyle=o_{\text{a.s. }}\left(n^{-r}\right),$ (52)
$\displaystyle\sup_{(u,v)\in[\delta,1-\delta]^{2}}\left|\bar{h}_{k_{e}\mid
j_{e};D_{e}}(u\mid v)-h_{k_{e}\mid j_{e};D_{e}}(u\mid v)\right|$
$\displaystyle=o_{\text{a.s. }}\left(n^{-r}\right).$

  * •

A3: For all $e\in E_{m},m=1,\ldots,d-1$, it holds:
(a) for all $(u,v)\in(0,1)^{2}$,

$\widehat{c}_{j_{e},k_{e};D_{e}}(u,v)-\bar{c}_{j_{e},k_{e};D_{e}}(u,v)=O_{p}\left(a_{e,n}\right),$
(53)

(b) for every $\delta\in(0,0.5]$,

$\displaystyle\sup_{(u,v)\in[\delta,1-\delta]^{2}}\left|\widehat{h}_{j_{e}\mid
k_{e};D_{e}}(u\mid v)-\bar{h}_{j_{e}\mid k_{e};D_{e}}(u\mid
v)\right|=O_{\text{a.s. }}\left(a_{e,n}\right),$ (54)
$\displaystyle\sup_{(u,v)\in[\delta,1-\delta]^{2}}\left|\widehat{h}_{k_{e}\mid
j_{e};D_{e}}(u\mid v)-\bar{h}_{k_{e}\mid j_{j};D_{e}}(u\mid
v)\right|=O_{\text{a.s. }}\left(a_{e,n}\right),$

where

$a_{e,n}:=\sup_{i=1,\ldots,n}\left|\widehat{U}_{j_{e}\mid
D_{e}}^{(i)}-U_{j_{e}\mid D_{e}}^{(i)}\right|+\left|\widehat{U}_{k_{e}\mid
D_{e}}^{(i)}-U_{k_{e}\mid D_{e}}^{(i)}\right|.$ (55)

  * •

A4: For all $e\in E_{m},m=1,\ldots,d-1$, the pair copula densities
${c}_{j_{e},k_{e};D_{e}}$ are continuously differentiable on $(0,1)^{2}$.

We note that A1 is satisfied by our marginal predictive estimator, as proved
in D.1 while A2, A3, A4 are all satisfied by the KDE pair copula estimator, as
shown in Nagler and Czado [2016].

#### Proof strategy

We perform the proof in three parts. To obtain the final result, we first
prove the convergence of pseudo observations to true observations through
induction. We then rely on the aforementioned convergence to show that
feasible pair-copula density estimators $\hat{c}_{j_{e},k_{e};D_{e}}$, and
conditional distribution function estimators $\hat{F}_{j_{e}|D_{e}}$ and
$\hat{F}_{k_{e}|D_{e}}$ are pointwise consistent. Lastly, these two results
are combined to obtain the convergence of the joint copula estimator
$\hat{\mathbf{c}}$ to the true multivariate copula $\mathbf{c}$.

#### Part 1: Convergence of pseudo observations

we start by proving a convergence rate of samples on the copula space obtained
through marginal distributions to their true unobserved equivalent. That is,
$\forall$ $e\in E_{1},\ldots,E_{d-1},i=1,\ldots,n$,

$\widehat{U}_{j_{e}\mid D_{e}}^{(i)}-U_{j_{e}\mid
D_{e}}^{(i)}=\mathcal{O}_{p}\left(n^{-r}\right),\quad\widehat{U}_{k_{e}\mid
D_{e}}^{(i)}-U_{k_{e}\mid D_{e}}^{(i)}=\mathcal{O}_{p}\left(n^{-r}\right).$
(56)

Starting with $e\in E_{1}$ (the conditioning set $D_{e}$ is empty), as a
consequence of A1 we obtain the bound

$\left|\widehat{U}_{j_{e}}^{(i)}-U_{j_{e}}^{(i)}\right|=\left|\widehat{F}\left(X_{j_{e}}\right)-F\left(X_{j_{e}}\right)\right|\leq\sup_{x_{j_{e}}\in\Omega_{X_{je}}}\left|\widehat{F}\left(x_{j_{e}}\right)-\overline{F}\left(x_{j_{e}}\right)\right|=\mathcal{O}_{p}\left(n^{-r}\right).$
(57)

To obtain the second part of (56) one can use identical arguments, providing
the initial inductive hook. Next, assuming (56) holds for all $e\in E_{m}$
with $1\leq m\leq d-2$, we extend the induction to $e\in E_{m+1}$. Recalling
that pseudo-observations $e^{\prime}\in E_{m+1}$ are equal to
$\hat{U}_{j_{e}|D_{e}\cup k_{e}}^{(i)}$ or $\hat{U}_{k_{e}|D_{e}\cup
j_{e}}^{(i)}$ for some $e\in E_{m}$, it follows by multiple triangle
inequalities that

$\displaystyle\bigl{|}\hat{U}_{j_{e}|D_{e}\cup k_{e}}^{(i)}-U_{j_{e}|D_{e}\cup
k_{e}}^{(i)}\bigr{|}$ $\displaystyle=$
$\displaystyle\bigl{|}\hat{h}_{j_{e}|k_{e};D_{e}}\bigl{\\{}\hat{U}_{j_{e}|D_{e}}^{(i)}|\hat{U}_{k_{e}|D_{e}}^{(i)}\bigr{\\}}-h_{j_{e}|k_{e};D_{e}}\bigl{\\{}U_{j_{e}|D_{e}}^{(i)}|U_{k_{e}|D_{e}}^{(i)}\bigr{)}\\}\bigr{|}$
$\displaystyle\leq$
$\displaystyle\phantom{+\;}\bigl{|}\hat{h}_{j_{e}|k_{e};D_{e}}\bigl{\\{}\hat{U}_{j_{e}|D_{e}}^{(i)}|\hat{U}_{k_{e}|D_{e}}^{(i)}\bigr{\\}}-\overline{h}_{j_{e}|k_{e};D_{e}}\bigl{\\{}\hat{U}_{j_{e}|D_{e}}^{(i)}|\hat{U}_{k_{e}|D_{e}}^{(i)}\bigr{\\}}\bigr{|}$
$\displaystyle+\;\bigl{|}\overline{h}_{j_{e}|k_{e};D_{e}}\bigl{\\{}\hat{U}_{j_{e}|D_{e}}^{(i)}|\hat{U}_{k_{e}|D_{e}}^{(i)}\bigr{\\}}-h_{j_{e}|k_{e};D_{e}}\bigl{\\{}\hat{U}_{j_{e}|D_{e}}^{(i)}|\hat{U}_{k_{e}|D_{e}}^{(i)}\bigr{\\}}\bigr{|}$
$\displaystyle+\;\bigl{|}h_{j_{e}|k_{e};D_{e}}\bigl{\\{}\hat{U}_{j_{e}|D_{e}}^{(i)}|\hat{U}_{k_{e}|D_{e}}^{(i)}\bigr{\\}}-h_{j_{e}|k_{e};D_{e}}\bigl{\\{}U_{j_{e}|D_{e}}^{(i)}|U_{k_{e}|D_{e}}^{(i)}\bigr{\\}}\bigr{|}$
$\displaystyle=$ $\displaystyle H_{1,n}+H_{2,n}+H_{3,n}$

Notice that for
$\delta_{i}:=\min\bigl{\\{}U_{j_{e}|D_{e}}^{(i)},U_{k_{e}|D_{e}}^{(i)},1-U_{j_{e}|D_{e}}^{(i)},1-U_{k_{e}|D_{e}}^{(i)}\bigr{\\}}>0$,
we have that all realisations $(U_{j_{e}|D_{e}}^{(i)},U_{k_{e}|D_{e}}^{(i)})$
are contained within $[\delta_{i},1-\delta_{i}]^{2}$ almost surely. Similarly,
all realisations $(\hat{U}_{j_{e}|D_{e}}^{(i)},\hat{U}_{k_{e}|D_{e}}^{(i)})$
are in $[\delta_{i}/2,1-\delta_{i}/2]^{2}$ for sufficiently large $n$ as a
consequence of (56). Combining this with A2 (b) and A3 (b), with large enough
$n$:

$\displaystyle H_{1,n}$
$\displaystyle\leq\sup_{(u,v)\in[\delta_{i}/2,1-\delta_{i}/2]^{2}}\bigl{|}\hat{h}_{j_{e}|k_{e};D_{e}}(u|v)-\overline{h}_{j_{e}|k_{e};D_{e}}(u|v)\bigr{|}=\mathcal{O}_{p}(a_{e,n}),$
$\displaystyle H_{2,n}$
$\displaystyle\leq\sup_{(u,v)\in[\delta_{i}/2,1-\delta_{i}/2]^{2}}\bigl{|}\overline{h}_{j_{e}|k_{e};D_{e}}(u|v)-h_{j_{e}|k_{e};D_{e}}(u|v)\bigr{|}=\mathcal{O}_{p}(n^{-r}),$

and by another application of (56),

$\displaystyle
a_{e,n}=\sup_{i=1,\dots,n}|\hat{U}_{j_{e}|D_{e}}^{(i)}-U_{j_{e}|D_{e}}^{(i)}\bigr{|}+\bigl{|}\hat{U}_{k_{e}|D_{e}}^{(i)}-U_{k_{e}|D_{e}}^{(i)}\bigr{|}=\mathcal{O}_{p}(n^{-r}),$

giving $H_{1,n}=\mathcal{O}_{p}(n^{-r})$. To complete part 1, we want to show
that $H_{3,n}=\mathcal{O}_{p}(n^{-r})$. We write the gradient of
$h_{j_{e}|k_{e};D_{e}}$ as $\nabla h_{j_{e}|k_{e};D_{e}}$ and use a first-
order Taylor approximation of
$h_{j_{e}|k_{e};D_{e}}\bigl{(}\hat{U}_{j_{e}|D_{e}}^{(i)}|\hat{U}_{k_{e}|D_{e}}^{(i)}\bigr{)}$
around $\bigl{(}U_{j_{e}|D_{e}}^{(i)},U_{k_{e}|D_{e}}^{(i)}\bigr{)}$ to get

$\displaystyle H_{3,n}$
$\displaystyle\leq\biggl{|}\nabla^{\top}h_{j_{e}|k_{e};D_{e}}\bigl{(}U_{j_{e}|D_{e}}^{(i)}|U_{k_{e}|D_{e}}^{(i)}\bigr{)}\begin{pmatrix}\hat{U}_{j_{e}|D_{e}}^{(i)}-U_{j_{e}|D_{e}}^{(i)}\\\
\hat{U}_{k_{e}|D_{e}}^{(i)}-U_{k_{e}|D_{e}}^{(i)}\end{pmatrix}\biggr{|}+o_{a.s.}\begin{pmatrix}\hat{U}_{j_{e}|D_{e}}^{(i)}-U_{j_{e}|D_{e}}^{(i)}\\\
\hat{U}_{k_{e}|D_{e}}^{(i)}-U_{k_{e}|D_{e}}^{(i)}\end{pmatrix}$

getting the desired result,. and hence the first equality of (56) by yet
another application of (56). The second equation follows by identical steps,
completing the induction.

#### Part 2: consistency of conditional CDF and pair-copula density estimators

Following similar steps to as in Part 1, one can obtain that for all $e\in
E_{1},\dots,E_{d-1}$, and all $\bm{x}\in\Omega_{\bm{X}}$, the CDF estimators
are bounded as

$\displaystyle\begin{aligned}
\hat{F}_{j_{e}|D_{e}}\bigl{(}x_{j_{e}}|{\bm{x}}_{D_{e}}\bigr{)}-F_{j_{e}|D_{e}}\bigl{(}x_{j_{e}}|\bm{x}_{D_{e}}\bigr{)}&=\mathcal{O}_{p}(n^{-r}),\\\
\hat{F}_{k_{e}|D_{e}}\bigl{(}x_{k_{e}}|{\bm{x}}_{D_{e}}\bigr{)}-F_{k_{e}|D_{e}}\bigl{(}x_{k_{e}}|\bm{x}_{D_{e}}\bigr{)}&=\mathcal{O}_{p}(n^{-r}).\end{aligned}$
(58)

To bound pair-copula density estimators, we apply the triangle inequality to
obtain

$\displaystyle\bigl{|}\hat{c}_{j_{e},k_{e};D_{e}}\bigl{(}u,v\bigr{)}-c_{j_{e},k_{e};D_{e}}\bigl{(}u,v\bigr{)}\bigr{|}$
(59) $\displaystyle\leq$
$\displaystyle\bigl{|}\hat{c}_{j_{e},k_{e};D_{e}}\bigl{(}u,v\bigr{)}-\overline{c}_{j_{e},k_{e};D_{e}}\bigl{(}u,v\bigr{)}\bigr{|}+\bigl{|}\overline{c}_{j_{e},k_{e};D_{e}}\bigl{(}u,v\bigr{)}-c_{j_{e},k_{e};D_{e}}\bigl{(}u,v\bigr{)}\bigr{|}$
(60) $\displaystyle=$ $\displaystyle R_{n,1}+R_{n,2}.$ (61)

Assumption A3 (a) coupled with (56) bounds $R_{n,1}$ while $R_{n,2}$ is
bounded byA2 (a), completing the second part.

#### Part 3: Consistency of the vine copula estimator

Up to now, our steps have mirrored those of Nagler and Czado [2016]. With the
following we differentiate ourselves by noticing that to get a bound on the
copula estimator alone, no marginal densities are required:

$\displaystyle\allowdisplaybreaks[1]\hat{\mathbf{c}}(\bm{x})$
$\displaystyle=\prod_{k=1}^{d-1}\prod_{e\in
E_{k}}\hat{c}_{j_{e},k_{e};D_{e}}\bigl{\\{}\hat{F}_{j_{e}|D_{e}}(x_{j_{e}}|{\bm{x}}_{D_{e}}),\,\hat{F}_{k_{e}|D_{e}}(x_{k_{e}}|{\bm{x}}_{D_{e}})\bigr{\\}}$
$\displaystyle=\prod_{k=1}^{d-1}\prod_{e\in
E_{k}}\biggl{[}c_{j_{e},k_{e};D_{e}}\bigl{\\{}\hat{F}_{j_{e}|D_{e}}(x_{j_{e}}|{\bm{x}}_{D_{e}}),\,\hat{F}_{k_{e}|D_{e}}(x_{k_{e}}|{\bm{x}}_{D_{e}})\bigr{\\}}+\mathcal{O}_{p}(n^{-r})\biggr{]}$
$\displaystyle=\prod_{k=1}^{d-1}\prod_{e\in
E_{k}}\biggl{[}c_{j_{e},k_{e};D_{e}}\bigl{\\{}F_{j_{e}|D_{e}}(x_{j_{e}}|{\bm{x}}_{D_{e}}),\,F_{k_{e}|D_{e}}(x_{k_{e}}|{\bm{x}}_{D_{e}})\bigr{\\}}+\mathcal{O}_{p}(n^{-r})+\mathcal{O}_{p}(n^{-r})\biggr{]}$
$\displaystyle=\mathbf{c}(\bm{x})+\mathcal{O}_{p}(n^{-r}).$ ∎

where the first line is a consequence of pair-copula estimator convergence in
Part 2, and the second equality is a consequence of (58) and the fact that
$c_{j_{e},k_{e};D_{e}}$ is continuously differentiable. This concludes the
proof.

##  Appendix E Experiments and practical details

Table 3: Average LPS (in bpd, lower is better) over five runs with standard errors for the Digits dataset.  Model | DIGITS
---|---
n/d | 1797/64
MAF | ${-8.76_{\pm 0.10}}$
RQ-NSF | ${-6.17_{\pm 0.13}}$
R-BP | ${-8.80_{\pm 0.00}}$
Rd-BP | ${-7.46_{\pm 0.12}}$
AR-BP | ${-8.66_{\pm 0.03}}$
ARd-BP | ${-7.46_{\pm 0.18}}$
ARnet-BP | ${-7.72_{\pm 0.28}}$
QB-Vine (30) | ${-6.56_{\pm 0.13}}$
QB-Vine (50) | ${-8.19_{\pm 0.14}}$
QB-Vine (100) | ${-8.54_{\pm 0.17}}$
QB-Vine (200) | ${-9.03_{\pm 0.17}}$
QB-Vine (300) | ${-10.24_{\pm 0.14}}$
QB-Vine (400) | ${-10.43_{\pm 0.13}}$
QB-Vine (500) | $\bm{-10.75_{\pm 0.03}}$
QB-Vine (full) | ${-10.11_{\pm 0.19}}$

Details of the UCI datasets are discussed in Ghalebikesabi et al. [2022]. In
experiments, We use the implementation of vine copulas from Schepsmeier et al.
[2015] through a Python interface. We follow the data pre-processing of
Ghalebikesabi et al. [2022] to make results comparable.

#### Hyperparameter search

In our experiments on small UCI datasets, we use a grid search over 50 values
from 0.1 to 0.99 to select $\rho$, independently across dimensions, selecting
possibly different values for each. To select the KDE pair copula bandwidth we
use a 10-fold cross-validation to evaluate the energy score for 50 values
between 2 and 4, as these ranges were appraised to give the best fits on
preliminary runs on train data. We note the hyperparameter selection of $\rho$
also supports gradient-based optimisation (see Appendix C). For energy score
evaluations, with marginal predictives we sample 100 observations and compare
them to the training data, while for the copula we simulate 100 samples from
the joint to compare with the energy score against training data.
For the PRticle filter we took an initial sample size of $d\cdot n$ to
accommodate for different dimensions while not being overcome by computational
burden. The Kernel used is a standard multivariate Gaussian kernel.

#### Compute

We ran all experiments on an Intel(R) Core(TM) i7-9700 Processor. In total our
experiments for the QB-Vine took a combined 15 hours with parallelisation
across 8 cores, or 120 hours on a single core. The Digits dataset on 8 cores
took us 6 hours to run with 5 different train and test splits. Other datasets
require about half an hour for five runs in parallel, while the Gaussian
Mixture study had a total time of 4 hours. The PRticle Filter takes about two
hours on all density estimation tasks combined. The RQ-NSF experiments on
Gaussian Mixture Models took about 4 hours combined. Our total compute time is
therefore the equivalent of 126 hours on a single core. Our implementation of
the QB-Vine is not fully efficient so the computational times are rough upper
bounds.

#### Selection of $P^{(0)}$ in practice

In practice, the initial choice of $P^{(0)}$ is made to reflect the support
and spread of the data. As we standardize our data to be mean $0$ and have
standard deviation $1$, a natural choice is the standard Gaussian
$\mathcal{N}(0,1)$. However, given distribution transformations are used
throughout the recursion, if observations fall in the tails of the predictive
density, numerical overflow might make them redundant, lowering the accuracy
of our approach. Therefore, it is desirable to have heavier tails than those
of the true distribution to capture outliers accurately. This coincides with
the theory on such recursion requiring heavier tails for the initial
predictive compared to those of the data, see the assumptions on $P^{(0)}$ in
Tokdar et al. [2009], Martin [2021] and Hahn et al. [2018], Fong et al.
[2023]. As such, our default choice is a standard Cauchy distribution.
For experiments, we tested Normal, Cauchy, and Uniform (over the range of the
training samples plus a margin) initial guesses. Generally, the Cauchy
distribution is a well-performing choice and obtained the best NLL in all but
two experiments. We give a summary of initial density $p_{0}$ choices for
different experiments in Table 4 for density estimation and Table 5 for
regression and classification.

Dataset | WINE | BREAST | PARKIN | IONO | BOSTON
---|---|---|---|---|---
Choice of $p_{0}$ | Cauchy | Cauchy | Cauchy | Normal | Cauchy
Table 4: Choice of $p_{0}$ for different density estimation experiments. Dataset | BOSTON (reg) | CONCR (reg) | DIAB (reg) | IONO (class) | PARKIN (class) |
---|---|---|---|---|---|---
Choice of $p_{0}$ | Normal | Cauchy | Cauchy | Cauchy | Cauchy |
Table 5: Choice of $p_{0}$ for different regression and classification
experiments.

#### Marginal Sampling

Here we briefly introduce the basic idea of the inverse sampling method for
marginal predictive densities via linear interpolation. In general, suppose
$Z\sim\mathbb{F}$ with the corresponding cumulative distribution function $F$,
we denote $F^{-1}$ as the generalized inverse function of $F$, given
$u\sim\mathcal{U}[0,1]$, then we can compute $z=F^{-1}(u)$ as a sample from
$\mathbb{F}$.
Suppose we have access to the cdf of a univariate random variable, meaning
that we can evaluate the cdf $F$ value given any location on the axis, and we
have a set of samples $z_{1:n}\stackrel{{\scriptstyle
i.i.d.}}{{\sim}}\mathbb{F}$. Our goal is to approximate the generalized
inverse $F^{-1}$ cdf then we can sample the realisations of this univariate
random variable as many as we want. Additionally, we assume that the random
variable is bounded and the closed support is
$\mathcal{R}=[\text{min}(z_{1:n})-\eta,\text{max}(z_{1:n})+\eta]$ where $\eta$
is a positive hyperparameter we manually tune as the range of extrapolation.
The idea is to construct a context set which contains $(K+2)$ pairs of
evaluation value on the axis and its associated cdf value and then we can use
this set to fit a linear interpolator to approximate the inverse cdf of our
desired distribution. More precisely, we can use a fixed grid of $K$ points,
$\bar{\bm{y}_{K}}=\\{\bar{y}_{1},...,\bar{y}_{K}\\}$ where
$\bar{y}_{1}<...<\bar{y}_{K}$, within $\mathcal{R}$ and compute
$\bar{c}_{i}=F(\bar{y}_{i})$ for $i=1,...,K$. Additionally, we set

$\begin{split}\bar{y}_{0}=\text{min}(z_{1:n})-\eta\,,\quad&\bar{c}_{0}=0\,,\\\
\bar{y}_{K+1}=\text{max}(z_{1:n})+\eta\,,\quad&\bar{c}_{K+1}=1\,.\end{split}$

Then, we obtain a context set
$\mathcal{C}(\mathcal{R})=\\{(\bar{y}_{i},\bar{c}_{i})\\}_{i=0}^{K+1}$. Since
the cdf is strictly non-decreasing and bounded within $[0,1]$, then we have
that

$\bar{c}_{0}\leq\bar{c}_{1}\leq...\leq\bar{c}_{K}\leq\bar{c}_{K+1}\,.$

Next, we propose to fit the approximation of the generalized inverse cdf
$\widehat{F}^{-1}$ via the linear interpolation using the context set
$\mathcal{C}(\mathcal{R})$ we constructed. Mathematically, given a value
$\bar{c}\in(0,1)$ we want to evaluate, then

$\bar{y}=\widehat{F}^{-1}(\bar{c})=\sum_{j=0}^{K}\frac{\bar{y}_{j}(\bar{c}_{j+1}-\bar{c})+\bar{y}_{j+1}(\bar{c}-\bar{c}_{j})}{\bar{c}_{j+1}-\bar{c}_{j}}\cdot\delta_{\\{\bar{c}_{j}<\bar{c}\leq\bar{c}_{j+1}\\}}\sim\mathbb{F}\,.$

By consequence, we can easily obtain the gradient of $\bar{y}$ w.r.t.
$\bar{c}$ as

$\frac{d\bar{y}}{d\bar{c}}=\sum_{j=0}^{K}\frac{\bar{y}_{j+1}-\bar{y}_{j}}{\bar{c}_{j+1}-\bar{c}_{j}}\cdot\delta_{\\{\bar{c}_{j}<\bar{c}\leq\bar{c}_{j+1}\\}}\,.$

Since we can evaluate the cdf values via $P_{j}^{(n)}$ for $j=1,...,d$, then
we can use the linear interpolation method we proposed. For
$\forall\,j\in{1,...,d}$, given the set of observations $x^{(1:n)}_{j}$, we
define the support for the $j^{\text{th}}$ marginal as

$\mathcal{R}_{j}=\Bigl{[}\text{min}(x^{(1:n)}_{j})-\eta,\text{max}(x^{(1:n)}_{j})+\eta\Bigr{]}\,.$

Following a similar procedure, we can construct a context set
$\mathcal{C}^{K}(\mathcal{R}_{j})$ with a grid of $K$ points via
$P_{j}^{(n)}$. Next, we can fit a linear interpolator denoted as
$\mathcal{LI}_{j}(\,\cdot\,)$ using the context set
$\mathcal{C}^{K}(\mathcal{R}_{j})$. Then,

$\mathcal{LI}_{j}(\epsilon)\sim
P_{j}^{(n)}\,,\quad\text{where}\;\epsilon\sim\mathcal{U}[0,1]\,.$

###  E.1 Comparison to normalising flow on Gaussian Mixture Model

We assess the performance of the QB-Vine on a mixture of 4 non-isotropic
Gaussians across a range of dimensions and sample sizes. We simulate $n$
$d$-dimensional data points from

$\displaystyle
p(\bm{y})\,=\,\sum_{k=1}^{4}\pi_{k}\cdot\phi(\bm{y};\bm{\mu}_{k},\bm{\Sigma}_{k})\,,$

where $(\pi_{1},\pi_{2},\pi_{3},\pi_{4})=(0.2,0.3,0.1,0.4)$ and

$\displaystyle\bm{\mu}_{k}\stackrel{{\scriptstyle
i.i.d.}}{{\sim}}\mathcal{U}[-50,50]^{d}\,,\quad\bm{\Sigma}_{k}\stackrel{{\scriptstyle
i.i.d.}}{{\sim}}\operatorname{Wishart}(d,\bm{I}_{d})\,.$

We compare the QB-Vine with the RQ-NSF as a benchmark off-the-shelf estimator.
The hyperparameters for the RQ-NSF were chosen to give the best performance on
train data, to be 100000 epochs, 0.0001 learning rate, 1 flow step, 8 bins, 2
blocks, and 0.2 dropout probability in common. For the number of hidden
features, we set $16$ for $d=10$, $32$ for $d=30$, $64$ for $d=50$, and $128$
for $d=100$. Our results in Table 6 show that the QB-Vine consistently
outperforms the RQ-NSF for the dimensions and sample sizes considered.

Table 6: Comparison of LPS for QB-Vine (our method) and RQ-NSF on GMM with 4 clusters for changing $n$ and $d$. Results for our QB-Vine method are shown as the top numbers of each row, and RQ-NSF values as the bottom numbers of each row. d \ n | 50 | 100 | 300 | 500 | $10^{3}$
---|---|---|---|---|---
10 | $\mathbf{{3.98_{\pm{0.23}}}}$ | $\mathbf{{1.73_{\pm{0.29}}}}$ | $\mathbf{2.15_{\pm{0.06}}}$ | $\mathbf{0.94_{\pm{0.31}}}$ | $\mathbf{2.43_{\pm{0.17}}}$
2-6 | ${36.47_{\pm{4.87}}}$ | ${17.14_{\pm{1.51}}}$ | ${12.82_{\pm{0.36}}}$ | ${7.10_{\pm{0.26}}}$ | ${7.91_{\pm{0.11}}}$
30 | - | $\mathbf{{17.94_{\pm{1.06}}}}$ | $\mathbf{{11.04_{\pm{0.35}}}}$ | $\mathbf{{12.87_{\pm{0.17}}}}$ | $\mathbf{{9.85_{\pm{0.40}}}}$
2-6 | - | ${91.09_{\pm{7.54}}}$ | ${50.51_{\pm{2.20}}}$ | ${48.50_{\pm{0.73}}}$ | ${34.98_{\pm{0.31}}}$
50 | - | - | $\mathbf{{38.59_{\pm{4.31}}}}$ | $\mathbf{{25.82_{\pm{0.06}}}}$ | $\mathbf{{26.14_{\pm{0.01}}}}$
2-6 | - | - | ${115.64_{\pm{3.06}}}$ | ${112.16_{\pm{2.05}}}$ | ${71.43_{\pm{1.65}}}$
100 | - | - | - | - | $\mathbf{{78.20_{\pm{0.23}}}}$
2-6 | - | - | - | - | ${268.88_{\pm{1.37}}}$