Quasi-Bayes meets Vines
Authors: David Huk
Yuanhe Zhang
Ritabrata Dutta
Mark Steel
Venue: NeurIPS 2024 poster (2024)
Keywords: Quasi-Bayesian
Copula
Vine Copula
Nonparametric Bayesian
density estimation
Primary area: probabilistic_methods
TL;DR: We introduced the Quasi-Bayesian Vine, which combines the strength of recursive joint Bayesian predictive densities with vines for varied tasks in high-dimensions, outperforming benchmark methods.

Abstract: Recently developed quasi-Bayesian (QB) methods \cite{fong2023martingale}  proposed a stimulating change of paradigm 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, however, existing constructions for higher dimensional densities are only possible by relying on restrictive assumptions on the model's multivariate structure. Here, we propose a wholly different approach 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. 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 \emph{an analytical form} and convergence rate independent of the dimension of the 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.

Quasi-Bayes meets Vines
David Huk
Yuanhe Zhang Department of Statistics
Department of Statistics University of Warwick
University of Warwick
David. Hukwwarwick.co ac.uk
Yuanhe. Zhang@varwick.ac.oo uk
Mark Steel
Ritabrata Dutta
Department of Statisties University of Warwick
Department of Statistics
University of Warwick
steel @warwick...... ac
Abstract
proposed stimulating Recently developed quasi-Bayesian (QB) methods
29
change of paradigm in Bayesian computation by directly constructing the Bayesian predictive distribution through recursion, removing the need for expensive compu- tations involved in sampling the Bayesian posterior distribution. This has proved be data-efficient for univariate predictions, however, existing constructions for higher dimensional densities are only possible by relying on restrictive assumptions on the model's multivariate structure. Here, we propose a wholly different ap- proach to extend Basss-Ba prediction to high dimensions through the use of Sklar's theorem, by decomposing the predictive distribution into one-dimensional predictive marginals and nigh-dimensiona copula. 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 fully non-parametric density estimator with an analytical form and convergence rate independent of the dimension of data some situations. Our experiments illustrate that the QB-Vine is appropriate for high dimensional distributions (~~4), needs very few samples to train (~2000) and outperforms state- of of-the-art art methods with analytical forms for density estimation and supervised tasks by a considcrable margin.
Introduction
The estimation of joint densities is a comeestonee of machine learning as a looking glass into the underlying data-generating process of multivariate data. Methods that support explicit density evaluation are crucial in probabilistic modelling, with applications in variational methods [58] 85).
Importance Sampling 2 (7), Sequential Monte Carlo [4,, Markov Chain Monte Carlo (MCMC)
92
49
and simulation-bait based inference [81 62 A prominent example Normalising Flows (NF)
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, with the Dirichlet Process Mixture Model (DPMM) [46] as the canonical nonparametric choice. Similar to kernel density estimation, the DPMM can be interpreted
80
24
infinite mixture of densities. It is composed of density called a kernel, with a fixed parametric form, and of a mixing density responsible for assigning those parameters. The specification of this kernel and mixing density is thus important
38th Conference on Neural Information Processing Systems (NcurIPS 2024).

---

as they regulate the expressivity of a DPMM. However, the computation of the DPMM's predictive density laborious relies on expensive MCMC methods scaling poorly high-dimension!!
Recently [9 proposed an appealing shift in Bayesian methods through the efficient constructions of Quasi-Bayesian (QB) sequences of predictive densities using only recursion, thereby circumventing the need for MCMC in Bayesian inference, and inspiring an offssring of works utilising this method-
ology 128
9 69 [OI]. The seminal univariate Recursive Bayesian Predictive (R-BP) its multivariate (R. extension [29] and the AutoRegressive Bayesian Predictive (AR-BP)
59 35
26
45
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.c. does not require large amounts of training data to be effective.
35
In multivariate QB predictives, to derive analytical expressions, it is necessary to enforce 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 similarity function of previous dimensions modelled through a covariance function or a deep autoregressive network, depending on the complexity of the data. The former relies on fixed form that is often too simplistic while the latter loses the appeal of data-efficien predictive like in the RB- Here, we posit that these assumptions, required for obtaining existing recursions, lack the flexibility to model multivariate data accurately.
In this paper, we introduce the multivariate QB recursions, utilising a copula decomposition to circumvent restrictive assumptions on the DPMM kernel. The QB-Vine is obtained by applying Sklar's theorem [9]] 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 simplified vine copula- highly flexible copula model suited for high dimensions. Compared to the sequential constructions of previous work, the QB- DB-e inherently parallelizable over dimensions instead of sequential. The main contributions of our work are as follows:
Vine (Qn-Vnne more general approach to forming
The copula decomposition frees us from the need to assume specific kernel structures of the DPMM, but preserves the data efficiency of QB methods while making more effective on high-dimensional data.
Under certain assumptions on the true dependence structure within the data, we show that the QB- attains a convergence rate that is independent of the dimension.
The above decomposition of the joint density and use of the energy score to tune hyperparam- eters makes the QB-Vine 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 [w we formulate the succinct survey of related work in Section and compare related methods to the QB-Vine in Section 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
Vine model. We provide
2 Quasi-Bayesian prediction
Notation. Let pxx) be a multivariate probability density function over X' C Rd, from which we observe i.i.d. samples Dp = {x ~ pxx). Similarly, let p² (x*)...... ......., (2d)) be the marginal densities of p(x), each over (a subset of) IR with corresponding ii.d. samples Dp. {z'x) p¹ (a*)) for i = 1,..., d, and assumed to all be continuous. Further, let P and pl pd be the respective cumulative distribution functions (cdfs) of the previously mentioned densities, in order of
For reference, with d the dimensionality of the parameter space. Random Walk Metropolis-HHatinne scales like OdP)) [87], Metropolis-adjusted Langevin algorithm like (dd/44) (86) and Hamiltonian Monte Carlo like (d433) M
2

---

appearance. Finally, when discussing predictive densities, we will use a subscript, plain for data (e.g. F,) and in parentheses for functions (e.g. P(n)) indicate the predictive at step n. distinguishing them from the superseript kept for dimension, and use bold fonts exclusively for multivariate objects.
Bayesian predictive densities as copula updates. Consider the univariate Bayesian predictive density P(n) for future observation R given seen i.i.d. data 21m E R" with likelihood of the data and a posterior ((n) for the model parameters 9 after n observations. By Bayes rule, we have:
((n)(T/Z1:)) = [ [(((0 f(xn 00) (mn (0 1im- d0
P(n 1) (2n (1:n-1)
As discovered by 43 29
50], multiplying and dividing by the predictive from the previous step
P(n-1j, we arrive
Joint density for x²
fxx(0) f(xn(0) T(n- (1 (f/x1:n-1) d8
P(n) = (
((n-1)(Znn (1n Pin 1) (z/(21n-1)
Marginal 2
Marginal
P(n (1 C(n)) (f'n-1) (x)), (xn))
which by Sklar's theorem 91 identifies the involvement of copulas in this recursive equation for the [ is seen to he a symmetric bivariate copula density function with the property that lim ((n) (2x, x²) = 1 a.s.Vr as a consequence of ((") with n [29]. It can be shown that every univariate Bayesian model can be written in this form 43), and so has unique copula sequence characterising its predictive updates by de Finetti's 's theorem 20 45 The choice of the predictive sequence then corresponds to
predictive densities. The second term on the right-hand side of
the almost sure convergence of p!"
an implicit choice of likelihood and prior 110 33
Marginal 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, [43] propose a nonparametric density estimator termed recursive Bayesian predictive (R-BP) Dirichlet Process Mixture Model (DPMM) inspired recursion emulating D They derive the correct Bayesian update under a DPMM for step and use it for all future steps m > (derivations shown in Appendix [B.3). The update copula of the R-BP is 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ío) and distribution P(o). the obtained analytical expression for the R-BP recursion has the following predictive density
(x) ()))) [(1
an) + an cp(P(n-1)xz) P(n-1)(2x.))
(2)
with c, E [0, 1] being bivariate Gaussian copula with covariance p. The corresponding edf also admits an analytical expression as follows
(r))()) (1 (a) (()-1((>) a. (,,,(((m-)))( P(n-)((n))
(3)
Here, for + the standard univariate Gaussian distribution,
Φ (u) p (c)
H,(,, u)
V1 p
is conditional Gaussian copula distribution with covariance p treated as hyperparameter and where Oe = (2 + J++T a sequence of weights converging to 0 (See supplement E of (29) for a more detailed explanation of the weights). The computational cost is O(n²) ) for initialising the recursion by computing P(R)(T))) for Ti:n and O(n) evaluate the pdf or edf point. This univariate R-BP model has been shown in [29 (Theorem 5) to converge in total variation to a limiting distribution P(oo) with density P(oo) termed the martingale posterior, defined in Appendix B which is quasi-Bayesiad object describing current uncertainty
29
31
50
8
33
43
shows that the
For an introduction to copulas, see Appendix
25 40

---

R-BP converges to the true density in the Kullback-Leibb distance. 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 8s algorithm, similar recursion for a DPMM's predictive [75] 744) ). Consequently, the R-BP 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 77 10 29 50) if the recursive updates are conditionally identically distributed (as is the case for the recursion of ()), they still exhibit desirable Bayesian characteristics such as coherence, regularization, and asymptotic exchangeability, motivating the Quasi-Bayesian Bayesian name as used in [30].
Multivariate recursive Bayesian predictives. An extension to the multivariate case was studied by [29] and refined by K5 for data with more complex dependence structures. The multivariate DPMM formulated as
/ K(x 0)dG(0), with G ~ DP Go) e
f(x
where K' is kernel for the observables x E E parameterised by 0, similarly to the kernel in kernel density estimation, and G is called the mixing distribution, upon which Dirichlet process prior is placed with base measure Gg and concentration parameter >> 0. To address this shortcoming, [29] provide a copula- based recursion obtained by assuming K ((x(0) = II (V 10", 1) and (0)
IV - 1), r > O, meaning both the kernel and base measure are assumed independent across dimensions, lacking the expressivity required to capture dependencies in the data. In [35], the
form of the kernel is relaxed be autoregressive with K (xx(0) = II=_, AV(x') 10"(1 (x) 1) where the kernel mean 0² : R + R is dependent on previous dimensions, and the base measure of [31] a product of Gaussian Process priors Go( = II GP(0 [0, 11k for k Ri-1 Ri-1 + R covariance function.
Vine copulas as effective models for high dimensions. Vine copulas are class of copulas that provide a divide-and-conquer approach to high-dimensional modelling by decomposing the joint copula density into d(d-1) bivariate copula terms. They considered among the best current copula models for density estimation. The main ingredient of a vine copula decomposition is the following identity consequence of Sklar's theorem 19]]:
p=l((=1aa") = c²b (P. Pb())) p=(x")
where a, b are subsets of dimensions from (1, (1.-..;00 d). Vine copulas rely on conditional factorisation p(x',...,
I 1 (x) to which they repeatedly apply 4 thereby splitting the joint density into the d marginal densities and d(d-1) bivariate copulas called pair copulas. The pair copulas for each i + j E {l.....dl .d., take as input pairs of conditional distributions (Pils* (zill*)) ) pulls (x=1s")) where S C (1....,d) d) (i,5) U is decided by the choice of the vine. A vine copula model thus has the form
II d(d-1)/2 d (pilsS (2415%) pils (191S") SS)).
(5)
iイj
We notice, that these pair copulas start as unconditional bivariate copulas and later capture higher orders of multivariate dependence by conditioning on the set S itself. This decomposition is valid but only unique up the permutation of indexes. We provide example of three-dimen dinessoam
vine copula decomposition and an overview in Appendix [.33, referring the reader to [17 [8] for
an introduction. In practice, we use simplified vine copula model 770] 71) which removes the conditional dependence of each of the copula e on SV, This is an approximation which reduces the complexity of the model for dependency structure but provides significant computational gains by reducing the size of the model space. An example simplified vine model in d = 3 is shown below: ((P)( (22)) ((3) ((())) = ::
- (P*(P((x))) (P"")))) (((3)) 22.1 (12 (18)(=)))))))))))
The number of pair copulas grows quadratically with the dimension, and the number of possible decompositions is exponential with the dimension, leading to greedy algorithms being used for model selection, see 22 The flexibility and efficacy of the vine has been studied in the literature, we refer to 17 18 73 93 among others. In particular, when the simplified vine assumption is true, [70] provides a dimension-inndpeedden convergence rate making simplified vine copulas greatly appealing for high-dimensional models.

---

3
Quasi-Bayesian Vine prediction
We propose the Buasi-Bayesian Vine (QB- (QB-Vine) for efficiently modelling high-dimensional pre- dictive density P(a) (and distribution P P(n)). ) The efficiency is achieved by adapting Sklar's theorem (19]). see Appendix A to split the joint predictive into predictive marginals for each dimension and high-dimensional copula:
Theorem 3.1. Let P(n) be an d-dimensional predictive distribution function with continuous marginal
(n)
predictive distributions P1 p2
pd (n) Then there exists copula distribution (n) such that
(m) (n)
for all x = (x", E ER::
x
P P(n)(x) ....., x²) ((.)(""....,="
(n)) (P))((())))))))) ..., Pa (m)((=))
And ifa probability density function available:
P(n)(22 ....,,"" x²)) = P(ng(a') (...)....... ... P(n)(x²) C(n) ((n) (Pmm)(")))........)) P(n)(
(6)
where ( ...) (1²) are the marginal predictive probability density functions, and C(n)
[0, 1] + R is the copula probability density function.
By applying the decomposition in (6) to two consecutive predictive densities P(m-1) and P(m)) we obtain a recursive update for joint predictive densities with two parts, of the form:
P(m) (x)
(am)( P Pmm--11 ()=))))))))))))))()))))
C(m) (Pimo (x ),.............).. P(i (m) (2²))
III
P P(m-1) (x)
C(m 1) (1 (.").... Pd (m- (x11)
Independent recursions
Implicit recursion on copulas
This decomposition fruitfully isolates updates to marginal predictive densities from updates their dependence structure, allowing us to model each recursion separately; the marginal predictives follow a univariate recursion a la while the copulas are free to follow another recursive form. Particularly:
As we are only interested in the joint predictive P(n). once marginal predictives are obtained, only need to fit single copula C(n) at step n to recover the joint predictive through (6)
where the recursion is done sequentially across dimensions, the QB-Vine can recurse marginal predictives in parallel by dimension.
Unlike 22 35
The model's dependence is not constrained by assumptions on the DPMM's form, instead, the QB- RB-Vine is free to fit any copula that best matches the dependence of the data.
Marginal predictive density estimation with the R-BP. We model marginal predictive densities P/n (Ir) and distributions Piny (Ix), vi 1,.... independently between dimensions. We the univariate R-BP approach described in Section 2 recursively obtain the analytical expression for both. For each dimension separately, starting with an initial density P/o; and distribution P we follow the updates for the density and (3) for the distribution.
Simplified vine copulas for high-dimensional dependence. After estimating marginal predictives, we model the joint density of (u² : p2 (2) (x)).......)))) == = pd (d (z with multivariate copula. We consider highly flexible simplified vine copula, found in Equation (5) which decomposes the
joint copula density c(u?,... ) into d(d-1) bivariate copulas " of the edfs from dimensions i and (possibly conditioned on additional dimensions) to capture the dependence structure of For the bivariate pair-copulas ctJ, we use a nonparametric Kernel Density Estimator (KDE)) Thus, each ed becomes two-dimensional KDE copula with the following expression:
Ei=
+ (u) --------iii (uL);
( Φ - (v)))
(u,v)
(x);0,b) (v-(())0,b)
where (b.; O, b) is the pdf of a normal with mean 0 and variance b, + is the inverse standard normal
edf. Samples [wu vi Jr_n casily obtained by iteratively fitting KDE pair copulas observed
samples {(u).-... .... u
For further detail KDE-based vine copulas, see Appendices A.2
A.3 and 70
5

---

Algorithm Joint, marginal, and copula density estimation with the Quasi-Bayesian Vine
Input: train M, N", d). teat(N, d), [BI...., Bi]. V.J.
Marginals:
1: for i = to do
# Dimensionnwwse training and evaluation, parallelised
2:
repeaa
3:
for m 1 to M do
I average over permutations, parallelised
4:
to N' Compute v.
5:
p (train") "P recursively
end for end for
7:
8:
Set P (N'+1> (-) s (N7 (1y (1e1 1/M
9:
Compute the energy score Sy Sp(1i. train') for min ~ P (1)()))))))))) using
10:
Update p² based Sg
11:
until convergenic of
12:
for to do
13:
Evaluate ((w141) ((train', ; teat'); recursively using
14:
Evaluate P( (ftrain', test") recursively using
15:
end for
16:
end for
Copula Fitting:
I: for b= B; Bi do
If Optimisation, parallelised
for == to N' mod(V) do
cross- valls-vaiiation, parallelised
3: Fit vine (N7 (a [
4: Compute the energy score (b(b) [u1. 1i (train lid
( ) for wiid vinel- using
5: end for
6: end for
7: Select 6" arg mins (D) for { [B.......
8: Evaluaae ((N'+a)(test) vinet Vinel (test) using
9: Evaluate P(yr- 1) (teat) = ViN' (test) IT Plwtt (teat'))
Return: P,Ny P1) P(N* (test), (N*+1) (test). pt.
Plad
Choice of Hyperparameters. We begin by choosing a Cauchy distribution the initial predictive distribution Pia) and the corresponding density P/o) Vi = 1....,d, d, with details provided in Appendix The hyperparameter p for the R-BP recursion is assumed different for each marginal predictive (i.e. p²)) Unlike previous R-BP works 43 29 351, the QB-Vine minimize the energy score to select p rather than the log-score, both of which are strictly proper scoring rules (see Appendix (aan (78) and define statistical divergences. This choice was motivated by the robustness properties of the energy score [77). The energy score is computed between observations and J 100 predictive samples from our marginal models conditional on previously observed data. As the R-BP sensitive to the ordering of the data, we follow 22 35] by averaging the resulting R-BP marginal over M 10 permutations of the data (see
for a discussion regarding the need of averaging over permutations).
We assume a same variance parameter for all the KDE pair copula estimators in the simplified vine and select it using 10-fold cross-validation, in data-dependent manner by minimizing the energy score between observations and J=100 copula samples. The assumption of a common bandwidth b motivated by mapping all pair copulas 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 (5). see Appendix (A.3) for which we use regular vine structure (79), selecting the best pair-copula decomposition with modified Bayesian Information Criterion suited for vines (72).
We include an algorithmic description for estimating the marginal density as well as the copula with the Qi-Vine in Algorithm where M is the number of permutations, N' and N the train and test sizes, d the dimension of the data, [B,..., Bi the copula variances considered for b, V the cross-validation size, and J the sample size for computing the energy score.
Computational benefits the Qi- Optimising the energy score instead of the log-score together with the copula decomposition, provides us with some significant computational gains.

---

Firstly, as the energy score is sample-based metric, we compute it through efficient univariate inverse probability sampling of P(m). Thus, we only recurse cdfs (()) during training, thereby halving the time to tune hyper-parameters compared to using the log-score which requires densities P(n) top of cdfs ((n) in (2). Secondly, Sklar's theorem implies the independence of marginal densities in the QB-Vine allowing to model them in parallel, reducing the cost to that of single R-BP, iee constant with d. Finally, the vine hyperparameter b is selected with a grid search using cross-validation, and parallelised across the grid and across cross-validation folds, each having the cost of a single vine.
Properties of the Quasi-Bayesinn Vine. To quantify the approximation of the QB- Vine, we provide the following stochastic boundedness [12] result for univariate R-BP distributions with respect to the limiting martingale posterior P(ax) (see [29] and Appendix B We note P(ae) is the univariate martingale posterior of the R-BP. The multivariate martingale posterior of the QB-Vine is not used in
our results, but we discuss properties, including a martingale condition, in Appendix [B.2
Lemma 3.2. (R-BP predictive distribution convergence) The error of the distribution function (x) in is stochastically bounded with
(P(s)(()) (in)(t)) = Op (1- wex
sup
Appendix D.I gives proof. In comparison to univariate KDE with mean-square optimal bandwidth 6, On (/5) which converges at rate Oa.s. (n 2/5 / //n(n)), the marginal R-BP has better rate with sample size. In what follows, we assume that the true copula is simplified vine of which
know the decomposition (a standard assumption in the vine copula literature [70] 18 (3)) 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, consider marginal distributions { { ((())) and {Pr (P))) are implicitly applied to x for respective copulas.
Theorem 3.3. (Convergence of Quasi-Bayesian Vine) Let (u) be the copula of {Pi (oo()))(")) )E=
and let C(n)(M) be the copula of (xin)(")){ Assuming that both copulas are simplified vine copulas with the same permutation of indexes in the decomposition of (). the estimation error is stochastically bounded V x E with
with
[((x)(x) (x)
where n is the convergence rate of the KDE pair-copula.
We provide proof in Appendix D.2 For bivariate KDE pair-copula estimator with optimal
bandwidth bn = O (n- 1//- we obtain n = n 1/3 (70]. From (95], we note the optimal convergence rate of nonparametric estimator is n (=/(/2++))) where is the number of times the estimator is differentiable. Therefore, as increases, we expect large benefits from using a vine copula decomposition for the QB-Vine. When the simplifying assumption does not hold, the simplified vine copula converges to a partial vine approximation of the true copula, as defined in [93]. Together, these two results guarantee accurate samples from the QB-Vine by inverse probability sampling arguments.
By Theorem 3.3 the copula C(m) ensures samples u = (u² (u),...,u ua) on the (0, 1j hypercuhe have dependence structure representative of the data. Then, marginal distributions Je recover dependent samples x E IRd by evaluating the inverse of the distribution at u² dimension-wise
Adapting the Vi-Vine for Regression/classifcatiio tasks. Our framework can accommodate regression and classification tasks in addition to density estimation by rewriting the conditional density as following:
xt p(x) = Pv(v)= 11 (p((2t) c(y, xt..,, (2 c(xl .d......))))
p(y|x) = ply. x)
cly,at,... PP(P))))) c(x² c(x)....,,2t)
The estimation of 7 is comprised of estimating the + marginals for y, xt rdd and the two copulas c(y.x x) and c(x). We specify separate P across marginal densities well separate
7

---

bandwidths b for the two copulas, to be estimated as in Algorithm I We note that for a copula decomposition be unique, we require that the marginals involved be continuous. This assumption is violated in the classification for binary outcomes As such, we make use of an approximation that transforms y to 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 [55] 61 44). 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 y = Ti+T where To and Ti are the numbers of negative and positive samples, respectively, in the training set. Consequently, we create different clusters in the copula [0, 1| 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 for classification with copulas which our framework can be extended to.
literature [1 15
18
4 Related Work
Our method shares similarities with existing work on QB predictive density estimation with analytical
forms. The pivotal works of [75 74] and the ensuing Predictive Recursion (PR) [37 65 66 97 63
36 64
propose recursive solution to the same problem but are restrained low dimensional settings due to the numerical integration of normalising constant over a space scaling with d. A sequential
importance sampling strategy for PR proposed 23 termed as PRticle Filter. The R-BP of (43
and the multivariate extensions in updates. In the multivariate case, imposing assumptions on the kernel structure leads to conditional factorisation of the joint predictive which recovers bivariate copula updates. In B5 an autoregressive Bayesian predictive (AR-BP) is used, where the dependence is captured by dimension-we -iiee 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 data-efficient predictive like in the R-BP. The Duasi-Bayesiaa 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-ding -igh-iimnnsional dependency through 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 [80 24 with analytical densities with a state-of-the-an performance across data types and tasks.
129 35 also have a recursive form driven by bivariate copula
5
Experiments
In this section, compare our QB-Vine model against competing methods supporting density evaluation with a bolosed-forma expression. Further details on the experiments are included in Appendix Code is included at htttps: http://gittuubt www/euk-oooooeeeeeo HooooHuuuDDvid/Du Vine
E
Density estimation. We evaluate the QB- on density estimation benchmark UCI datasets 4 with small sample sizes ranging from 89 to 506 and dimensionality varying from 12 to 30, adding results for the QB- and PRtcle Filter to the experiments of 35). We report the log predictive score LPS Riestt 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 (83), DPMM (84) with a diagonal (Diag) and full (Full) covariance matrix for each mixture component, MAF 80), RQ-NSF 24 as well as the closely related PRticle Filter 23), R-BP 29 and AR-BP B5 For the last two Bayesian predictive models, we add a subscript d to indicate that the o hyperparameter possibly differs across dimensions, and the net suffix indicates network-basse selection of p for dimensions. We observe in Table I that our QB-Vine method comfortably outperforms all competitors as the dimension increases, while getting close to the performance of the best alternative Bayesian predictive models for the lower dimensional WINE dataset. Our method's 's relative performance increases with the dimension of the data, particularly achieving much smaller LPS for IONO the dataset with the largest dimensions and relatively small sample size. We accredit this performance to the copula decomposition as that is our main distinguishing factor from the other Bayesian predictive models.
Σ Eesr' - II (Xx)) on a held-out test dataset of size Nteatt
8

---

Table I: Average log predictive score (lower is better) with error bars corresponding to standard deviations over five runs for density estimation on datasets analysed by 33 dimension increases, the Vine outperforms all benchmarks.
We note that as
WINE 89/12
BREAST 97/14
PARKIN 97/16
IONO 175/30
BOSTON 506/13
nld
KDE
13.69
10.45
12.83
32.06
8.34 0.00
DPMM (Diag)
17.46+0.00
16.26+0.71
22 22.2240.000
35.30+1.28
7.64+00.0
DPMM (Full)
32.88+0.82
26.67
39.95+1.56
86.18+10.22
9.45:20.43
MAF
39.60+1.41
10.13440.0
11.76+0.45
140.09
56.01 4:2.774
RQ-NSF
38.34
+6.57
31.26 31.264+0.51
54.49
-2 -2.20 22.00=0.11
PRicle Filter
23.89
25.9811.00
34.79
79.22219.66
27.18-33.12
R BP
13.57
7.45.0.02
9.15g0.04
21.15
4.5620.04
RR-RP
13.32+0.01
6.12=0.00
7.52=0.05
19.82=0.08
-13.50+0.50
AR-BP
13.45+0.00
6.18+0.05
8.29-0.11
17.16+0.26
-0.45-0.77
ARg BP
13.2210.04
6.11200.0
2.21
16.48-10.26
-14.7520.80
ARnet-BP
14.41+0.11
6.87 6..2
8.29+0.17
16.32+0..56
-5.71-0.62
QB Vine
13.76
4..7
0.
16.08 12
-31.0411.00 1.00
High- -iggeniiinni image dataset. We fur- ther evaluate the QB- QB-Vine on the digits dataset (n 1197, d=64) as high-dimensional exam- ple with a relatively low sample size. The high contrast between n and d makes the problem suited for assessing the data efficiency and con- vergence of the QB- Vine. We compare with the two NF models as their high model capacity is good fit for image data, as well all the Bayesian predictive methods of Section [5] from the study of [35]. We report the average LPS in bits per dimension (bpd) with standard errors over fifteen runs with random partitions, using half the sample size to train models and the other half evaluate the LPS. Additionally, we report the average LPS of the QB-Vine obtained in the same way except for the training set size be- ing reduced (to 30, 50, 100, 200, 300, 400, 500). F depicts the QB- Ri-ev performance for different-sized training sets. When trained on the full training 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 training size of 200, demonstrating its data-efficiency and convergence speed. A complete numerical table is reported in Appendix
QBVine
-7
RQ NSF
LPS in bpd
RD
AR-BP
-9
ARnet-UP
200 400
600
800
Training sample s22e
Figure 1: 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.
Regression and classification. We further demonstrate our method's effectiveness on super- vised learning tasks, with three datasets for regression and two datasets for classification, adding to the study of (35) For classification, we transform the binary values to continuous ones to preserve copula assumptions, as detailed in Section
3 We report the conditional LPS =
Wiest Ex=r - In (P( (irraat.))) (yk (x+)) over test set of size Thtas made up of half the samples with the conditional estimator trained on the other half of the data. We compare our model against a Gaus- sian Process [100). a linear Bayesian model (Linear) (68), one-hidden-lay multilayer perceptron (MLP), as well the R-BP and AR-BP variants for supervised tasks (29 35). The QB-Vine outperforms competing methods on all datasets except CONCR. We believe the lower performance on CONCR 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 ML.)), which we posit is the consequence of our conditional estimator in S defined as a ratio, inflating the variation in the LPS. However, we highlight that an overly precise inference is more misleading/dangerous than overly uncertain one.

---

Table 2: Average LPS (lower is better) with error bars corresponding to standard deviations over five runs for supervised tasks analysed by [5 The QB-Vine Vine performs favourably against benchmarks, with relative performance improving samples per dimension decrease.
Regression CONCR
Classification
BOSTON 506/13
DIAB
TONO
PARKIN
nid
1,0308
442/10
351/33
195/22
Linear
0.87++.08
0.99
1.07+0.01
0.33++..01
0.384
GP
0.42-4o.06
0.36 02
1.062...0
0.30
0.42:g0.02
MLP
1.4211.01
2,01 a8
3.32
0.261.0.0
0.31g0
R. BP
0.7610.00
0.87 0a
1.05
0.26
0.37
0.00+
0,42++0.0 00
1.00 02
0.34+0.02
0.27++0.03
AR BP
0.52.+o.13
0.42++0.01
1.06:=0.02
0.21+0.02
0.29::0.02
1.
0.39
0.99 0.00
0.20
0.28+0.03
ARnet-BP
0.45+0..11
0.03--0.00
1.41+0.07
0.24++004
0,26+0.04
Qn-Vine
0.81
0.54 84
0.87
-1.88=11.1
-0.76 28
Scalability of the Vine: In Appendix F.T we assess the scalability of the QB-Vine on large data sizes and dimensions, by fitting Gaussian mixture models with 20 random means and non- isotropic covariance in dimensions = 50 to d = 600 with n = 10000 train and test sets. We compare our model to RQ-NSF, reporting the LPS as well as the maximum mean discrepancy and the reverse Kullback-Leible divergence to assess sample quality, showing superior performance.
Discussion
We introduced the Quasi-Bayesian Vine, joint Bayesian predictive density estimator with an analytical form and easy to sample from. This extends the existing works on Quasi-Bayesian predictive densities, by using Sklar's theorem decompose the predictive density into predictive marginals and copula to model the high- dimenniienii dependency. This decomposition enables two-part estimation procedure, employing Baass-Bayesiann recursive density estimation for the marginals and fitting simplified vine copula for the dependence, resulting in convergence rate independent of dimension for certain joint densities. We empirically demonstrate the advantage of QB- Vine on 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 main bottleneck of the QB-Vine is the simplified vine, both computationally and methodologically. The non-uniqueness of the vine decomposition resulting in a search over an exponentially large model space during estimation and the hyperparameter selection of the KDE pair copulas could both lead misspecified models. Further, our main assumption is the use of simplified vine copula which is only an approximation to the true distribution. While these concerns stem from limited effective copula models for high dimensions being available, from practical point of view, simplified vine offers tractable and fast likelihood evaluations, and ultimately outperforms competitors as shown in experiments.
Future directions this work include the incorporation of more effective copula models, or copulas accommodating different dependence structures [98 73 51). Another exciting direction are develop- ments of new recursive quasi-Bayes methods that can be merged into a Quasi-Bayesing Vine model
331
Author Contributions: David Huk wrote the code, ran the experiments, derived the proofs and wrote the paper. Yuanhe Zhang helped write the initial code and the appendix. Ritabrata Dutta conceptualised the project and David Huk formulated the concrete method. Ritabrata Dutta and Mark Steel jointly supervised the project and helped to write the paper.
10

---

Ackmowledgments: We thank all the reviewers for their helpful feedback. We also thank Surya Tokdar, Fabrizio Leisen and Edwin Fong for useful discussions. further thanking Vaidehi Dixit for sharing Prticle Filter codes with us. David Huk is funded by the Center for Doctoral Training in Mathematical Sciences at Warwick. Ritabrata Dutta is funded by EPSRC (grant nos. EP/V025899/1 and EP/T017112/1 and NERC (grant no. NE/T00973X/1).
References
[1] Pierre Alquier, Badr-Eddine Chérief-Abdelaaiff, Alexis Derumigny, and Jean-David Fermanian. Estimation of copulas via maximum mean discrepancy. Journal of the American Statistical Association, pages 1-16, 2022.
[2] Michael Arbel, Alex Matthews, and Arnaud Doucet. Annealed flow transport Monte Carlo. In International Conference on Machine Learning, pages 318-330. PMLR, 2021.
[3] 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
/
//openrvview.coo wwt htte///ooounni-=ttdddk.....
[4] Arthur Asuncion and David Newman. UCI machine learning repository, 2007.
[5] Tim Bedford and Roger M Cooke. Probability density decomposition for conditionally depen- dent random variables modeled by vines. Annals of Mathematics and Artificial intelligence, 32:245-268, 2001.
[6] 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: arXiv:1705.10744 2017.
[7] Patrizia Berti, Luca Pratelli, and Pietro Rigo. Limit theorems for a class of identically distributed random variables. 2004.
[8] Patrizia Berti, Luca Pratelli, and Pietro Rigo. Almost sure weak convergence of random probability measures. Stochastics and Stochastics Reports, 78(2):91 2006.
[9] Patrizia Berti, Luca Pratelli, and Pietro Rigo. Limit theorems for empirical processes based on dependent data. 2012.
[10] Patrizia Berti, Emanuela Dreassi, Fabrizio Leisen, Luca Pratelli, and Pietro Rigo. A proba- bilistic view on predictive constructions for Bayesian learning. Statistical Science, 1(1):1-15, 2023.
[11] Alexandros Beskos, Natesh Pillai, Gareth Roberts, Jesus-Maria Sanz-Sema, and Andrew Stuart. Optimal tuning of the hybrid Monte Carlo algorithm. 2013.
[12] Yvonne Bishop, Stephen Fienberg, and Paul W Holland. Discrete multivariate analysis: Theory and practice. Springer Science & Business Media, 2007.
[13] Arthur Charpentier, Jean-David Fermanian, and Olivier Scaillet. The estimation of copulas: Theory and practice. Copulas: From theory to application in finance, 35, 2007.
[14] Yuhui Chen, A copula-based supervised learning classification for continuous and discrete data. Journal of Data Science, 14(4):769-782, 2016.
[15] Yuhui Chen and Timothy Hanson. Copula regression models for discrete and mixed bivariate responses. Journal of Statistical Theory and Practice, 11:515-530, 2017.
[16] David Clayton. 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.
[17] Claudia Czado. Analyzing dependent data with vine copulas. Lecture Notes in Statistics, Springer, 222, 2019.

---

[18] Claudia Czado and Thomas Nagler. Vine copula based modeling, Annual Review of Statisties and Its Application, 9:453-477, 2022.
[19] A Philip Dawid, Monica Musio, and Laura Ventura. Minimum scoring rule inference. Scandi- navian Journal of Statistics, 4331):123-138, 2016.
[20] Bruno De Finetti. La prévision: ses lois logiques, ses sources subjectives. In Annales de I' I'institut Henri Poincaré, volume 7, pages 1-68, 1937.
[21] Jeffrey Dissmann, Eike C Brechmann, Claudia Czado, and Dorota Kurowicka. Selecting and estimating regular vine copulae and application to financial returns. Computational Statistics & Data Analysis, 59:52-69, 2013.
[22] Vaidehi Dixit and Ryan Martin, Permutation- hased uncertainty quantification about a mixing distribution. arXiv preprint arXiv: XXXv::996.0034// 2019.
[23] Vaidehi Dixit and Ryan Martin. A PRticle filter algorithm for nonparametrie estimation of multivariate mixing distributions. Statistics and Computing. 33(4):1-14. 2023.
[24] Conor Durkan, Artur Bekasov, lain Murray, and George Papamakarios. Neural spline flows. Advances in neural information processing systems, 32, 2019.
[25] 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.
[26] Fabian Falck, Ziyu Wang, and Christopher C Holmes. Are large language models bayesian? martingale perspective in-context learning. In ICLR 2024 Workshop Trustworthy Large Language Models, 2024.
Secure and
[27] Matthias Fischer, Christian Köck, Stephan Schlüter, and Florian Weigert. An empirical analysis of multivariate copula models. Quantitative Finance, 9(7):839-88) 854, 2009.
[28] Edwin Fong and Brieuc Lehmann. predictive approach to bayesian nonparametric survival analysis. In International Conference on Artificial Intelligence and Statistics, pages 6990-7013. PMLR, 2022.
[29] Edwin Fong, Chris Holmes, and Stephen G Walker. Martingale posterior distributions. Journal of the Royal Statistical Society Series B: Statistical Methodology, 85(5):1357-1333 2023.
[30] 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-1144 2020.
[31] Sandra Fortini and Sonia Petrone. Prediction-based uncertainty quantification for exchangeable sequences. Philosophical Transactions of the Royal Society A, 381((227):202211422 2023.
[32] David Frazier, Jeremias Knoblauch, and Christopher Drovandi. The impact of loss estimation Gibbs measures. arXiv preprint arXh::2404.15649 2024.
[33] Samuele Garelli, Fabrizio Leisen, Luca Pratelli, and Pietro Rigo. Asymptotics of predictive distributions driven by sample means and variances. arXiv preprint arXiv:2403. 16828, 2024.
[34] Gery Geenens, Arthur Charpentier, and Davy Paindaveine. Probit transformation for nonpara- metric kernel estimation of the copula density. 2017.
[35] Sahra Ghalebikesabi, Chris C Holmes, Edwin Fong, and Brieuc Lehmann. nonparametric density estimation via autoregressive predictive updates. In Uncertainty in Artificial Intelligence, pages 658 -68. PMLR, 2023.
[36] Subhashis Ghosal and Aad Van der Vaart. Fundamentals of nonparametric Bayesian inference, volume 44. Cambridge University Press, 2017.
12

---

[37] Jayanta K Ghosh and Surya T Tokdda.. Convergence and consistency of Newton's 's algorithm for estimating mixing distribution. In Frontiers in statisties, pages 429 443. World Scientific, 2006.
[38] Tilmann Gneiting and Adrian E Raftery. Strictly proper scoring rules, prediction, and estima- tion. Journal of the American statistical Association, 102(477):359-3177 0785999-7777- 2007.
[39] 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.
[40] Joshua Größer and Ostap Okhrin. Copulae: An overview and recent developments. Wiley Interdisciplinary Reviews: Computational Statistics, e4(3) 2022.
[41] Shixiang Shane Gu, Zoubin Ghahramani, and Richard Turner. Neural adaptive sequential Monte Carlo. Advances in neural information processing systems, 28, 2015.
[42] Ingrid Hobæk Haff, Kjersti Aas, and Arnoldo Frigessi. On the simplified pair-copula construc- tion-simply useful or too simplistie? Journal of Multivariate Analysis, 10165):129- 2010.
[43] P Richard Hahn, Ryan Martin, and Stephen G Walker. On recursive Bayesian predictive distributions. Journal of the American Statistical Association, 1136523):1855-109; 2018.
[44] 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):e20022SS003 120, 2022.
[45] Edwin Hewitt and Leonard J Savage. Symmetric measures on Cartesian products. Transactions of the American Mathematical Society, 80(2):470-501, 1955.
[46] Nils Lid Hjort, Chris Holmes, Peter Müller, and Stephen G Walker. Bayesian nonparametrics, volume 28. Cambridge University Press, 2010.
[47] 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 2013.
[48] Marius Hofert, Ivan Kojadinovic, Martin Mächler, and Jun Yan. Elements of copula modeling with R. Springer, 2018
[49] Matthew Hoffman, Pavel Sountsov, Joshua V Dillon, Ian Langmore, Dustin Tran, and Srinivas Vasudevan. Neutralizing bad geometry Hamiltonian Monte Carlo using neural transport. arXiv preprint arXiv: 03704, 2019.
[50] Chris Holmes and Stephen G Walker. Statistical inference with exchangeability and martingales. Philosophical Transactions of the Royal Society A, 38 381(2247):2022144,, 2023.
[51] David Huk, Rilwan A Adewoyin, and Ritabrata Dutta. Probabilistic rainfall downscaling: Joint generalized neural models with censored spatial Gaussian copula. arXiv preprint arXiv:2308.0983 2023.
[52] David Huk, Lorenzo Pacchiardi, Ritabrata Dutta, and Mark Steel. David Huk, Lorenzo Pacchiardi, Ritabrata Dutta and Mark Steel's contribution the discussion of 'Martingale posterior distributions' by Fong, Holmes and Walker. Journal of the Royal Statistical Society Series B: Statistical Methodology, 85(5): 405-1406, 2023.
[53] Tim Janke, Mohamed Ghanmi, and Florian Steinke. Implicit generative copulas. Advances in Neural Information Processing Systems, 34:26028-26033 2021.
[54] Harry Joc. Dependence modeling with copulas. CRC press, 2014.
[55] André G Journel and Wenlong Xu. Posterior identification of histograms conditional to local data. Mathematical Geology, 26:323-359, 1994.
13

---

[56] Sanket Kamthe, Samuel Assefa, and Mare Deisenroth. Copula flows for synthetic data generation, arXiv preprint arXiv:2001.00598 2021.
[57] Leonid V Kantorovich. Mathematical methods of organizing and planning production. Man- agement science, 6(4):366-422) 1960.
[58] Diederik Kingma and Max Welling. Auto-encoding variational Bayes. arXiv preprint arXiv:/312.6114. 2013.
[59] Hyungi Lee, Eunggu Yun, Giung Nam, Edwin Fong, and Juho Lee. Martingale posterior neural processes. The Eleventh International Conference on Learning Representations.
[60] Chun Kai Ling, Fei Fang, and Zico Kolter. Deep Archimedean copulas. Advances in Neural Information Processing Systems, 33:1535-1545, 2020.
[61] Weiwei Liu. Copula multi-label learning. Advances in Neural Information Processing Systems, 32, 2019.
[62] Jan-Matthis Lueckmann, Jan Boelts, David Greenberg, Pedro Goncalves, and Jakob Macke. Benchmarking Aimulation-basee inference. In International conference artificial intelli- gence and statistics, pages 343-351, PMLR, 2021.
[63] Ryan Martin. Convergence rate for predictive recursion estimation of finite mixtures. Statistics
& Probability Letters, 82(2):378-335 384, 2012.
[64] Ryan Martin. A survey of nonparametric mixing density estimation via the predictive recursion algorithm. Sankhya B, 83:97-121, 2021.
[65] Ryan Martin and Jayanta K Ghosh. Stochastic approximation and Newton's s estimate of mixing distribution. Statistical Science, pages 365-382, 2008.
[66] Ryan Martin and Surya T Tokdar. Asymptotic properties of predictive recursion: robustness and rate of convergence. 2009.
[67] Alex Matthews, Michael Arbel, Danilo Jimenez. Rezende, and Arnaud Doucet. Continual repeated annealed flow transport Monte Carlo. In International Conference Machine Learning, pages 15196-15219, PMLR, 2022.
[68] Thomas Minka. Bayesian linear regression. Technical report, Citeseer, 2000.
[69] Blake Moya and Stephen G Walker. Full uncertainty analysis for bayesian nonparametric mixture models. Computational Statistics & Data Analysis, 189:107838, 2024.
[70] 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.
[71] Thomas Nagler, Christian Schellhase, and Claudia Czado. Nonparametric estimation of simplified vine copula models: comparison methods. Dependence Modeling, 5(1):99-120, 2017.
[72] 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.
[73] Thomas Nagler, Daniel Krüger, and Aleksey Min. Stationary vine copula models for multi- variate time series. Journal of Econometrics, 227(4)305-324,, 2022.
[74] Michael A Newton. On nonparametric recursive estimator of the mixing distribution. Sankhya: The Indian Journal of Statistics, Series A, pages 306-322, 2002.
[75] Michael A Newton, Fernando A Quintana, and Yunlei Zhang. Nonparametric Bayes methods using predictive updating. In Practical nonparametric and semiparumetric Bayesian statistics, pages 45-61 Springer, 1998.
14

---

[76] 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.
[77] Lorenzo Pacchiardi, Sherman Khoo, and Ritabrata Dutta. Generalized Bayesian likelihood-free inference. arXiv preprint arXiv:2104.0333 03889, 2021.
[78] Lorenzo Pacchiardi, Rilwan 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.
[79] Anastasios Panagiotelis, Claudia Czado, Harry Joe, and Jakob Stöber. Model selection for discrete regular vine copulas. Computational Statistics & Data Analysis, 106: 136-122,-, 2017.
[80] George Papamakarios, Theo Pavlakou, and lain Murray. Masked autoregressive flow for density estimation. Advances in neural information processing systems, 30, 2017.
[81] George Papamakarios, David Sterratt, and lain 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.
[82] 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.
[83] Emanuel Parzen. On estimation probability density function and mode. The annals of muthematical statistics, 33(3): 1962.
[84] Carl Rasmussen. The infinite Gaussian mixture model. Advances in Neural Information Processing Systems, 12, 1999.
[85] 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.
[86] Gareth O Roberts and Jeffrey S Rosenthal. Optimal scaling of discrete approximations to Langevin diffusions. Journal of the Royal Statistical Society: Series B (Statistical Methodol- ogy), 60(1))255-268, 1998.
[87] GO Roberts, Gelman, and WR Gilks. Weak convergence and optimal scaling of random walk Metropolis algorithms. The annals of applied probability, 7(1):110-122) 1997.
[88] 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 puckage version, 2(5), 2015.
[89] Jayaram Sethuraman. A constructive definition of Dirichlet priors. Statistica sinica, pages 639-650, 1994.
[90] 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.
[91] Sklar. Fonctions de répartition à dimensions et leurs marges. Annales de I'ISUP, volume 8, pages 229-231, 1959.
[92] Jiaming Song, Shengjia Zhao, and Stefano Ermon. A-nice-me: Adversarial training for meme. Advances in neural information processing systems, 30, 2017.
[93] Fabian Spanhel and Malte S Kurz. Simplified vine copula models: Approximations based on the simplifying assumption. Electronic Journal of Statistics, 13:1254-1291, 2019.
[94] Bharath K Sriperumbudur, Arthur Gretton, Kenji Fukumizu, Gert RG Lanckriet, and Bern- hard Schölkopf. A note on integral probability metrics and g-divergences. arXiv preprint arXiv:0901.2698, 2009.
15

---

[95] Charles Stone. Optimal rates of convergence for nonparametric estimators. The annals of Statistics, pages 1348-1360, 1980.
[96] Gábor J Székely and Maria L Rizzo. new test for multivariate normality, Journal of Multivariate Analysis, 93(1):58-80, 2005.
[97] Surya T Tokdar, Ryan Martin, and Jayanta K Ghosh. Consistency of recursive estimate of mixing distributions. The Annals of Statistics, pages 2502-2522, 2009.
[98] Thibault Vatter and Thomas Nagler. Generalized additive models for pair-copula constructions. Journal of Computational and Graphical Statistics, 27(4):715-727, 2018
[99] Ziyu Wang and Christopher C Holmes. On uncertainty quantification for near-bayes optimal algorithms. In Sixth Symposium on Advances in Approximate Bayesian Inference-Non Archival Track.
[100] Christopher Williams and Carl Rasmussen. Gaussian processes regression. Advances in neural information processing systems, 8, 1995.
[101] Luhuan Wu and Sinead A Williamson. Posterior uncertainty quantification in neural networks using data augmentation. International Conference on Artificial Intelligence and Statistics, pages 3376-3384. PMLR, 2024.
16

---

A Copulas
Copulas are a widely adopted tool in statistics and machine learning for modelling densities, permit- ting the construction of joint densities in 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.
Definition 2.1.1 in [8]] (Copula) copula distribution a multivariate distribution function with standard uniform univariate marginal distributions, that is, (O0,1) 1) margins.
This approach motivated through Sklar's theorem:
Theorem A.1 (Sklar (9I). Let P be d-dimensional distribution, function with continuous marginal distributions Pl ......, (x
,, Pd. Then there exists a copula distribution such that for all x
x)) E Rd:
P(x', ....., x²)) C(P) (. P,
P(a',...,
And if probability density function (pdf) is available:
p(a' ((.",..,<²) p² (x²) (x')).....'' (2*) c(Pl (x*)).... pa(
where p (x*),.... .e.)),..,,""" ,, are the marginal pdf', and (( (x*),.... pd P' the copula pdf.
If the marginal distributions are absolutely continuous, the copula is unique. Consequently, one can decompose the estimation problem learning p(x) by first learning all the marginals {or 1d Mez1. and in second step learning an appropriate copula model c(u', ..... ) where u' == Pi(x"), ) (1.... ... are the images of the x' under the edf of each dimension. By applying edf transformations marginally,
copula agnostic to the differences between dimensions such as axis scaling, and purely focuses on 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 [27 47). The Gaussian copula a popular parametric choice as it is well-studied and can be fitted quickly even to moderate dimensions [48]. However, it lacks the desired flexibility to capture more complex dependencies involving multiple dimensions. Among nonparametric copulas, Kernel Density Estimator (KDE) copulas [3 are commonly used. They apply an inverse Gaussian distribution the observed {u 1 map them to latent space and perform regular KDE on the latent density. However, this KDE copula method suffers from the poor scaling of KDE estimators 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 candidate solutions such as [60] [56] 53 3 As such, current copula models are mostly limited to low to
modelling (5]
A.1 Gaussian copula
popular parametric copula model is the Gaussian copula. assumes that the dependence between dimensions is identical to that of a Gaussian distribution with mean 0 and covariance matrix Σ;
c(u' wd) Na (u')..... 2- E)
IT N(o 0,
As such, its only parameters are the off diagonal entries of the covariance matrix E. In the case of = 2, there is a single parameter to estimate for the Gaussian copula.
A.2 Gaussian Kernel Density Estimator copulas
As the marginal distribution of copula is uniform in [0. 1]. the support of the copula estimator is restricted to J0, 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 13 34 is popular approach for such nonparametric copula models, which models the copula on latent Gaussian space. We explain the approach the following.
17

---

By the inverse sampling theorem, for any R-valued continuous random variable X, applying the corresponding cumulative distribution function F to X results in F(X) being uniformly distributed in [o, 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 the inverse standard normal distribution, then the latent distribution for each marginal will be normal and 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 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:
c(u, v) E @ (w)
+ + T (ux); o - (v=)).0))
お Φ (u):0,)
o @ (o)(0))
where (u, v) and (ux, Ub) are both in [0, 1|2 [ug, 2k k== observed copula samples, and & and s are respectively the Gaussian distribution and density with mean 0 and variance h > 0.
A.3 Vine copulas
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 [5]. Here we illustrate the decomposition by a 3-dimensional copula density cu,vw for the random vector (U, V, W):
(u, = cu. crw(y,,)) co,wiv (Cuyv(u) w), Cw (ucl (w)lw)
where
Cu.v is the copula of (U, V), Cv.w is the copula of (V, W);
Cuwivv is the copula of (Colv (U]u), Cwiv (Wl v)) conditional on V == u;
Criy is the conditional distribution of Cu,v on V, and Cw is the conditional distribution of Cy,w on V.
Generally, the distribution of CU,WiV= 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:
cu.,ww U, w) = cu.v (u, v) cv, ww) cu.iiv (Culv(ul Cwiv(wl
The rationale of the simplified assumption is studied in 42). In this paper, we mainly focus on the regular vine copula (R-vine) with simplified assumption. The construction an R-vine copula has two basic ingredients: (1) valid tree structure of all random variables, (2) the choice of family for bivariate copulas.
Before we introduce the tree structure which valid construct R-vine copula, let us first rigorously define the random variable of a copula we want to estimate. Suppose U= (U1. U2. Ua) is d-dimensional random variable which is distributed as a copula C. then we have that each marginal random variable U; is uniformly distributed in [0. 1]-scale. For notational simplicity, we will use the index of random variable as its notation instead. Denote T = (Ti,... Ta-1) } as sequence of (d - 1) trees in terms of T (Vi, &.). Here we use set of two nodes to represent the corresponding edge in T., i.e., e = (a, b) node and in ε, are linked. To construct valid R-vine copula, T satisfies the following conditions:
18

---

Ti is tree with a set of edges E₁ and set of nodes V1 = (I,..., d]. dt:
T is tree with set of edges & and set of nodes V; = -(E;-1) for 2,,,..., 1, where VEE) denoted that the pair of nodes which are linked by an edge in E is treated new node.
If two nodes in Est1 are linked by an edge in V4+1, then they must be linked to one common node in Ti.
For Ve = () E To with 2 2, we define
de) = a
b(e) - b(
Finally, we rigorously define the R-vine copula follows.
Definition A.2 (Regular Vine Copulas). A d-dimensional copula C' is regular vine copula if there exists a uple (T,C) such that
T is regular vine tree sequence with (d 1) trees;
C {Ca E E [d 1], C. is bivariate copula} is family of bivariate copulas for each edge;
For Ve E &, and Vi E [d - 1], then C. is corresponding to the copula of (ate). c(e)
Therefore, the density function of R-vine can be expressed as
C(T,C) (M1, (u1, Ud)
d-1
II cEE.
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.
CT ((1, U2, U3 u4, us) = =c(11, U2) c(u1, uz) c(u2, u4) c(us, us)
11,5 [1j3 Us}s) C2,3/ i (U2 1 Uzn) (1,4/2(1122 uaj2)
2,,5 (22)1,3 U5J1,3) 03,4/1,2 2,,41,2(23331122 (#3)1,2, 21122)
2,,,51,,,3 (44 1,2,3 U511,2,)
Specifying the tree structure for R-vine decomposition is essential and plays an important role in pair-copula estimation through conditioning. An excellent overview is given in [17]. Notably, an R-vine decomposition is not unique for given joint copula pdf. An appealing tree selection algorithm is proposed in (72), where authors derive modified BIC criterion that prioritizes sparse trees while being consistent when the dimension d grows most at the rate of yn where is the sample size.
B Martingale Posterior Distributions
Here we explain martingale posterior distributions as justification Bayesian approach through focus on prediction.
B.1 The Bayesian Choice as Consequence of Predictive Uncertainty
common goal statisties is the inference a parameter or quantity 0 by analysing data in the form of observations (f1...., Ix), E N. The rationale for learning from data is that each observation provides information about underlying process and parameter without which statistical analysis would be redundant. Indeed, consider decision-maker basing their decision their belief about (in an i.i.d. setting). Having observed data (æ1,...,, 2nn and given the opportunity to observe an additional data point Zn+1, they would be assumed to accept, they could refine their beliefs
19

---

on A. 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 9 captured by the prior (0). the goal is the inference about the distribution of the parameter given observed data (z1,..., ,nn) which is encoded in the posterior density [(n) (f((2),..., T..)
For the decision-makes to refuse, it would mean that the additional observation has no significant effect on their belief. This implies, as identified by (2) and [50], that there is point where additional observations (In+1,... IN) provide no benefit to the knowledge update of 0. Inspecting the
Bayesian posterior in terms of observed data 21n = (a1....., and possible future observations to be made 2n+1: written as
((n) (0) f ([n)) (0, 100 d1n++1000
one can expand the right-hand integrand by including the predictive density for future data, obtaining
T(n) (0) (01) = f (cc) (0)21:0)) predictive density
1:00+1100/210 1:00 (21:))
Bayes
Having rewritten the posterior density in this way, it becomes apparent that the uncertainty in the value of f, which is given by the Bayes estimate, is a consequence of the uncertainty surrounding the imputation of missing observations In+1: through the predictive. With this insight. unlike the traditional Bayesian construction of p with predictive mechanism using all available data to impute the missing 1n+1: (or at least impute FnN for a sufficiently large N), and replace the Bayes estimate π ((x)) with an appropriate functional of the complete data Z1:00- This predictive mechanism used impute the unobserved data
[29] proceed to replace the predictive density
1:01:00 given observed I.m directly leads to the martingale posterior as the limiting distribution of functionals when the unobserved data has been imputed, see definition (29..
B.2 The martingale posterior of the QB-Vine
the authors are interested in predictive resampling (not considered in this paper for the QB-eine), which means to progressively sample from P(n) and using those samples instead of observations to continue the recursive construction of P(++1)- For this predictive resampling to make sense (see conditions and in (2)), they show that it is sufficient for the sequence to be conditionally identically distributed. This condition is in turn shown to be equivalent to the following martingale condition in [9:
In (29 35|
Theorem B.1. (Theorem 3.1 of (9]). A sequence of random variables X1, X2,.. is c.i.d. and only if its densities P2....... are such that for n 2 0 and all x E Rd:
/ P(n)(x) P(n-1)(x")dxx -P(n-1)xx).
We show that our QB-Vine construction indeed satisfies this condition. The only assumption the ratio of two consecutive vines must be I, i.e., the dependence structure is constant between predictive steps. We note that these derivations hold for any marginal recursive construction of the form (2) and any copula density used for x,, but show it for the QB- -i-Vine here. We interpret the condition of having the same dependence structure between steps as natural when data is supposed to come from the same data-generating process, which is indeed the circumstance in which we apply the QB-Vine. Further, given observations, the best guess of the multivariate copula is given by fitting it at the last iteration, which our approach in practice. We write c(0) = v() to avoid confusion with bivariate copulas and denote inputs to v()) as vectors (Pi (xr1)... Pa(aa)) = (Pt())) For
clarity, we also write the mixture between the independence and Gaussian copula used in (
as [(1 - an)) + On cpp (P (2)) P(n-1) (x.))) = (P(2211) (x*), P( (a()). )). The proof as
follows:
20

---

Proof.
/ ((x)(x)) P(==1)(2")dd""
vo (PP(((( P(n- P(n-1)(xx) dx"
d
(P(o-1)(x') ((an)((Pp-))))))))
(m)((Pim) (x')))
II cd C(P Pl P(m-))(()))))) 1=1
v(a-1)( ((--)(((n-_)(=))]] ((n-)((P(nn))) (لقيع) dx
1
(P v(a)(PPi)) (x'))) ((--1)([u)) du"
-P(n-1)(x)
V(n (111)(11(00 (1)(()))))(()))
-Pin-1)(X)
V(n-1) (111
(z"), x/)) ((ux)) v6) - a,) -(nm-))(t)) a (((()))()))) du"
V(n) (m(((inn)(( -P(m-1)(X) V(n-1) (1--))(P(n (1P((11))) P(n-1)(xX.
The first equality is applying Sklar on (x) and using recursion on the ensuing marginal densities (a²). The second equality is Sklar on P(n-1)(Xx.) The third step is obtained by the
substitution duf, Pin) (#t,)dz Lastly, we use equation (3) for the edf inside the copula. Then,
the result follows by noticing that the bivariate copulas and the copula ((n=1)([(,,)) integrate to
1 by copula properties (see e.g. the proof of Theorem 6 in 9) and [(1 -a,) P (1
(ar)
H (PPin) (u,)] integrates into [P] (x="))("") due to it being a martingale marginally for each i. O
Consequently, predictive resampling is also possible with our approach and is as simple as sampling from the fitted copula, and marginally updating each univariate R-BP. This can be done in parallel across dimensions, instead of sequentially as in [29 33), which is computationally appealing and opens an interesting avenue for future work.
B.3 Constructing martingale posterior for the DPMM
As explained in Section 2 given sequence of observations (2:)121 from the data generating process, the Bayesian approach wants obtain an infinite sequence of copula densities (G))21 recursively update the Bayesian predictive distribution. However, finding such an infinite sequence of copulas isn't 't feasible in practice and the ensuing copula family is determined by the likelliod pair one implicitly chooses. For example, from [3). we have the following copula sequences:
If we select the likelihood and prior as ((y; 0) = De By and ग (0) = e -- then (n 1) u)
n (1 -1 1/n (1 vn)) 1/1 -
Cn (u,u, vn)
1|
which sequence of Clayton copula [16] with parameter
If we select the likelihood and prior as ()y )) = o(y; 0 1) and T()) = 600;0,7 then
On (u, vn) = 2₂ ( (u), 4 (v,) Op.))
o (u): 1) 0,
where Φ2 is the pdf of a bivariate normal distribution, Φ is the pdf of normal distribution, and
1
Σ,n
Pn. = (n 71
Pn
21

---

Therefore, we obtain a sequence of Gaussian copulas with parameters (p h221. Further, if we additionally put conjugate prior on the variance parameter o² of the likelihood 1(v(0) = (vu; 0, r²), then we will recover sequence of Student-t copulas.
The issue of model mis-specification naturally arises here due to the selection of prior-likelihodd pair. Fortunately, here we can employ the DPMM which has relatively high flexibilityy due to its nonparametric nature,For the DPMM, suppose the first observation ø1 arrives, then we can obtain the first predictive via the updating kernel
E (f) Sa(a1) [[f((=)) [[fo(=1)]
kn (1) =
where po( B[fa(x)]. Then, we can derive the numerator as E [fa(2) fc(#1)
= E
Wy we o 0z; 1) (x1; 0x. 1)
E
Wi
Eca [0(x; 0, 1)] EG. (0(z1 0,
+
EG. jo(z; 0, 1) ((11; 0,
Σ
po(a) po( +E
Ec, [0(x; 0, 1) c(x); 0, 1)]
The first equality follows from the stick-breaking representation
89 of the DP, we can formulate G
as
G( =
Ws 50, (-)
where
w. = U; II 0z),
2/2 3.dd Beta(1,
0,
fid. G
Then, the second equality follows from the condition that Σ D Wi almost surely. Denote EE
|Σ w then we can write the updating kernel
(x, = (1 a) a EGo (6(:; 0, 1) d(a1;00001 1)] Po(0) Po(
22 4 (c),
= (1 a) +
σ (w):0, 6 -v1);0,,
(1 a)
spl $ (w), Φ (v1)),
where
1 covariance matrix Σ,- P
1
T- is the inverse edf of standard normal distribution,
Po(x), U1 Po(()))
po( = -000
22

---

is a bivariate Gaussian copula density function with the correlation parameter 11 T).
Therefore,
that the copula ei (u, 01) for the first updating step is mixture of independent copula and Gaussian copula. However, we will lose the tractable form of the copula from the second updating step. According to [43], instead of deriving the explicit rule for the sequence copula, we fix the correlation parameter and set a to be a (0, 1)-valued decreasing sequence (a.)221
C Strictly Proper Scoring Rules
In probabilistic machine learning. a Scoring Rule (SR) measures the appropriateness distribution in modelling observation E X' through score written S(P. x). For (hyper) parameter estimation, suppose we aim to model the underlying distribution an observation x using a family of distributions I' parametrized by 0, then we can SR S to select the suitable 0. If we assume that x ~ Q, then we can obtain the expected SR S via taking expectation w.r.t x as
S(P, Q
Ex~o (s(p
According to 38 we call S strictly proper S(P, Q) is minimized if and only if = Q. i.e. for VP P with Q such that
S(Q. Q) < S(P, Q)
Ltd. Considering our settings of marginal distributions, we have observed data X ~ IP" and we aim to model IP using Pa. More explicitly, we don't need to use Pd directly in the expected SR, instead we normally use its probability density function or samples from this distribution to evaluate. In general, optimise
8"
ag Q)
Age
argmin S(P., Q)
AEe
Since we do not have the complete population of Q we use the empirical SR S instead, i.e.
O = S Q)
868
where S(Po, Q) E S(F. Under mild conditions, it can be proven that 0- + 0-
asymptotically in 19
For any positive definite kernel k( -, -), the Kernel Score (38 is given by
(T(Pox)) E Y') -2 (Y,(
where Y, Y Pp. If we set k(x, y)
- |x yllg- then we obtain the energy score which is strictly proper Ey-P, Y < co. The energy Score defined as
S (Po. = 2 E Y
E Y Y for B E (0,2).
Here YandY' are i.i.d. samples from Pa. Practically, given finite samples y (0), ..., ym (0) i.d Pa where y0) denotes sample from Pe which is differentiable function of 0, the unbiased estimate of the energy Score via Monte Carlo approximation is
3 (y (0),
(0)
x 8
(1y, (0)
(8)
m(m 1)
Similarly, we can derive the unbiased gradient of Se (y1.m (0), x) w.r.t. which is crucial for any gradient descent algorithm. For notational convenience, define
(Y,Y = x Y Y
23

---

then S (Po.x) P Y', ,x)]. Next, we have that
VV Y ( Y, Y xx]
=.y Vog Y. Y' x)
(0),),(0)))xx0))) x) 61331)
m(m
2 (0) x11 m(m -1 k#j
Volly, yx(0)1
m
S (y1.m) 1:m (0), x)
Furthermore, if P and Q both are univariate distributions of R-valued random variables, then the energy Score will reduced to
E (P) = 2 E - x/ - - Y for 3 E (0,
(9)
where Y,Y" ~ Po.. where the expectations can similarly be approximated using samples from Fa. Notice that this will become the Continuous Ranked Probability Score (CRPS) 196] when B = 1.
The energy Score is strictly proper scoring rule [38 19 76 and is special instance of the maximum mean discrepancy 39] as well as a statistical divergence. has been used as an effective objective for copula estimation 53 51 and even for R-BP marginal predictives
shown in [52]. It has also enjoyed success as an objective for Normalising Flows [90] and generative models 778 The energy Score is the only objective among Wasserstein p-metrics [57] that supports unbiased gradient evaluations [6, (6). has known optimisation-free solution and features faster convergence rate than similar integral probability metries (94 32
D Proofs
D.1 Lemma 3.2
Proof. We follow the Definition 14.4-3 of stochastic boundedness from [2]. Begin by choosing
8 E (0, 1). From Proposition [29]. we have for and M > n, over the supremum of E X::
2
P ((P((a)) P(m) 2 e) S 2exp
20a
Σ
Di=n+1
a
2
lim P ((o)a( se)
lim 1- 2exp M-000
M-+o
2ea
+
1 Σ Ei=n+1
a
Next, we choose a value e (where the subscript shows that this quantity is dependent on n) to have the appropriate probability the right-hand side by enforcing:
2
8 = 2 exp
26+0 3
a
i=n+1
= En + log
8
2a
En +
1
log
Σ 1=++1
00
a
3
2
with solution
20.+1 12
- log (9) 20
log
2 2log ($)
De Eient1 oz
en
Due to 8 (0, 1) and a; Ovi we have:
log (3) 12/22 2an+1 2 3
log (9) n 1/2 20
+
- 2 log (3) E a?
3
[, 112
24

---

With the choice of a; (2
) we see that limn- n aa is bounded for all powers
2 -. As such, both log (2) n 1/2222n+ and log (3) n n/12 20 will safely be bounded
for large enough n. Similarly, due to the choice of di we have Li (a,) O(n meaning n n Eien 2i=+t(aa) (ax)² will be bounded for large enough n long as -1. Hence -2 log (5) n Lin++1 of will also be bounded for large enough n..
Consequently, for our choice of 8, there exists a finite K > 0 and a finite N > 0 such that:
sup P (1 (Mioole) P(n) (I) =EX'
Vn > N. S 2 1 1 δ
n n-1/2 1/2
D.2 Theorem
33
We prove the statement for general densities, which then naturally extends predictive densities. The proof of this result is largely an adaptation of the simplified vine copula convergence result (Theorem in (7]) but where no rate on marginal densities required. As such, our proof shares an identical approach until the last part, where we deviate. We have 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 superscripts with parenthesis to now differentiate between samples instead of predictive steps, following notational conventions of the literature. We define l--functions the conditional distribution functions for pair copulas
Ie (u) = 1 Cjele::: (s, v)ds, for (u, v) E [0, 1).
Further, we refer to the unobserved samples of pair-copulas as
(rl joID. Ffe De
De (x([) B.).
(-(i] DD, FR.ID. (x
D.
(10)
for = 1,..., n. We also denote estimators and quantities obtained by application to these unobserved samples with bar superscript, for example:
De (u, u) Cic (u, (1) ... j.......... U..ID.) (...)))))............))))
Finally, denote with hat superscript all quantities and estimators obtained by using 0j (0) instead of the true unobserved samples used in Equation (TO)
())
Assumptions For completeness, we state assumptions about the marginal distribution estimator P well bivariate copula estimators è, even though our practical choices from the main text respect these. We begin by stating the assumption about our marginal distribution estimator denoted P:
Al: The marginal distribution function estimator has the following convergence rate:
sup P(x) - P(z) Op n 112 zEXt
Next, we state our assumptions about the pair copula estimator for completeness:
A2: For all E Fm m I....., 1, with -r the convergence rate of a bivariate KDE copula estimator, it holds: (a) for all (u, v) 1).
De..D. (u, (u, Ce.kee (u, (n-"),
(b) for every 8 (0,0.5).
h h Mj.lko (u v) - hj (u u) Oas. (n-))
sup (u,))
I (u v) - h Dalli (u u) Oas. (n (r-))
sup
(6.1))
25

---

A3: For all e E Em, m = ......d. . it holds:
(a) for all (u, v) 1)²
Decke: Decke:De (u, Op (dc,))
(b) for every 8 (0,0.5).
sup (u.+)e(6.1-8f
Io.DD (u (0) - hs v)
Oas. (aem)
T. (u v) hin (u u) Oas. (Dem)
sup
(v.v)e(6.1-683
where
((0)
De - ((()
| 1 | - -U(()
sup
J.ID.
KeID.
i=....,7
A4: For all e Em: differentiable on (0, 1)
d the pair copula densities DD., continuously
We note that A1 is satisfied by our marginal predictive estimator, as proved in D.I while A2, A3, A4 all satisfied by the KDE pair copula estimator, shown in (70]
Proof strategy We perform the proof 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 &.D... and conditional distribution function estimators F2. D. and Fs.IDD. are pointwise consistent. Lastly, these results are combined obtain the convergence of the joint copula estimator the true multivariate copula c.
Part 1: Convergence of pseudo observations we start by proving convergence rate of samples on the copula space obtained through marginal distributions to their true unobserved equivalent. That is, V Ei,..., E.-1. 1,...,7. n,
(
((()) DD =O, (n
000 U D.|De D.
Op (n-").
(11)
D.
Starting with E Ej (the conditioning set Dr is empty), as consequence of A1 we obtain the bound
10 v = i -F
F
sup zy E0x
(F = Op (n
To obtain the second part of (1] one can use identical arguments, providing the initial inductive hook. Next, assuming 11 holds for all e E Em, with 1 S m S d - 2, we extend the induction to E Em+1. Recalling that pseudo-observations e' Em+1 are equal to ((()) some e Em, it follows by multiple triangle inequalities that
(+())
for j.lD.uk. k.IDD.jj.
10:00
J.,,,Uk - - D.ukll
[haakee: ( (6)
j.ID, K.ID. hilBB.D. (0 10 (0
10,0
j.lo.
k.I...
[ (0,00
D. (0)
D. Tylk..D. D. (0,0) joD,
(0) (000 K.
Mu. ( (6(1)) (0()) D.D -..lk... (0(0) 1000
10(0
D. ko D.
[hiekee { (6(1)
D. (0(0) p. h (T(() 10()
JalD. [
H1, + H2n + Hs,n
u(() tions (U(1) (7(1) DeID. ke D. j=|D. ke D. > 0, we have that all realisa-
Notice that for 5, == min (1))
ul)
(1)
D k.) are contained within 55, 5,12 almost surely. Similarly, all realisations
D.
26

---

(0 J.I.. k.I.D. are in
for sufficiently large consequence of (1
Combining
this with A2 (b) and A3 (b), with large enough
1 (2.0)6(6,/211-6.33))
sup
(t.l::D.(ul))) (ulu) = Op(den):
H2.n (m.)))(64/21116//2
sup
[njolks:D. (u/e) - bg = Op(n
and by another application of
0g = sup 10:00 U(1) 10 i=1..... D. jolD, K.ID. K.ID. 0n
U
giving Hinn Op(n-)). To complete part 1, we want to show that H3.n - Op(n~)) We write the gradient of MelkeDD as Vhjalke:D. and use a first-order Taylor approximation of ((() 1.1 k.) around (U ja|D.> koID.) to get
Di,,,, (Ú(() ((
ja|D.
(i()) (=(()
(t
[((1
H3, & VTT.IIR..DD (U (
De D.)
Jel
D.
+ Oa.s.
D. D2
De
(rl1) U
((0)
U(()
keID.
kelD.
Ro
Re
D.
getting the desired result,, and hence the first equality of (1I) by yet another application of [1] The second equation follows by identical steps, completing the induction.
Part 2: consistency of conditional CDF and pair-copula density estimators Following similar
steps as in Part 1, one can obtain that for all EI...., Ed-1, and all estimators are bounded as
Qxx the CDF
P.. (.D)) Fj.I (xD.) = Op(n
(12)
F..ID. (Iv, (TD.) - FK.ID. (Ix, (.p.) = Op(n To bound pair-copula density estimators, apply the triangle inequality obtain
1eje.re:D. (u, Decke v))
< (ej (u,v) De (u, + I... (u, CjakeiDe (u,
RA,1 Rn,2-
Assumption A3 (a) coupled with bounds Rr,1 while Rn is bounded byA2 (a), completing the second part.
Part 3: Consistency of the vine copula estimator Up to now, steps have mirrored those of 770. With the following we differentiate ourselves by noticing that to get a bound on the copula estimator alone, marginal densities are required:
d-1
c(x)
II A=1 E
II
Ce. {Ei D. (.),).DD)) (o,
II E
Se.haiD (xD) Pr. D. (ox, (D.)) (p(n)r)
II
II
Sj.DD.DD (F (25, (DD.), Fk. D. (ex, (D.)) Op(n Op(n
A= EE
c(x) + O,,nn
where the first line is a consequence of pair-copula estimator convergence in Part 2, and the second equality is a consequence of (12) and the fact that &.... is continuously differentiable. This concludes the proof.
27

---

Table 3: Average LPS (in bpd, lower is better) over five runs with standard errors for the Digits dataset.
Model n/d
DIGITS 1797/64
MAF
-..6.10.10
RQ-NSF
-6.17=0.13
R-BP
-8.80.0.00
R.-BP
-7.46=0.12
AR-BP
-8.66
ARg-BP
07.46
ARnet-BP
- -7.72=0.28
QB- (30)
6.66
QB- QB-Vine (50)
-7.92=0.22
QB-Vine (100)
-8.98
QB- Vine (200)
-9.69=0.48
QB- (300)
-10.23
QB-Vine (400)
-10.39+0.13
QB-Vine (500) -10.1940.20
QB-Vine (full)
-10.47
E Experiments and practical details
Details of the UCI datasets are discussed [5 In experiments. We use the implementation of vine
copulas from [8] through Python interface. We follow the data pre-pro pre-processing of 33 to make results comparable.
Hyperparameter search In experiments small UCI datasets, values from 0.1 to 0.99 to select p. independently across dimensions, selecting possibly different values for each. To select the KDE pair copula bandwidth we use 10-fold validaaiit 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 p also supports gradient-based optimisation (see Appendix C which we utilise in the Gaussian Mixture Model general, gradient-baa
grid search 50
optimisation of p converges within less than five iterations. For energy score evaluations, with marginal predictives, we sample 100 observations and compare them to training data, while for the copula we simulate 100 samples from the joint to compare with the energy score against training data.
experiments Appendix[E.]
For the PRticle filter, we took an initial sample size of d accommodate for different dimensions while not being overcome by computational burden. The Kernel used is standard multivariate Gaussian kernel.
Compute We ran all experiments ments for the QB-Vine took a combined 15 hours with parallelisation across cores, or 120 hours on single core. The Digits dataset 8 cores took 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 total time of 4 hours. The PRticle Filter takes about two hours all density estimation tasks combined. The RQ-NSF experiments on Gaussian Mixture Models took about hours combined. Our total compute time is therefore the equivalent of 126 hours single core. Our implementation of the Di-e is not fully efficient so the computational times are rough upper bounds.
Intel(R) Core(TM) i7-9700 Processor. In total our experi-
Selection of P(o in practice In practice, the initial choice of P(o) is made 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 (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 (()) in 99 64 and 43 a standard Cauchy distribution.
29
As such, our default choice
28

---

For experiments, we tested Normal, Cauchy, and Uniform (over the range of the training samples plus margin) initial guesses on train data. Generally, the Cauchy distribution is werf-oornonnnng choice and obtained the best NLL in all but two experiments. We give a summary of initial density po choices for different experiments in Table for density estimation and Table 5 for regression and classification.
Dataset
WINE
BREAST
PARKIN
IONO
BOSTON
Choice of Po
Cauchy
Cauchy
Cauchy
Normal
Cauchy
Table Choice of po for different density estimation experiments.
Dataset
BOSTON
CONCR (reg)
DIAB
IONO (class)
PARKIN (class)
Choice po
Normal
Cauchy
Cauchty
Cauchy
Cauchy
Table 5: Choice of po for different regression and classification experiments.
Marginal Sampling Here we briefly introduce our inverse sampling method for marginal predictive distributions via linear interpolation. In general, through basic rules of probability, for univariate æ ~ p with the corresponding cumulative distribution function P with inverse p-1, given (000, 1]. we can obtain ā = P=-(u) as a sample from p. However, in our case, we have an analytical expression for P only, with no expression for p-1 As such, to sample from our model, need approximation of p-1 that we can evaluate.
To do so, we start by considering the range of values over which we seek to approximate P-1, In our work, we consider a range = [min( - v) max(#1:.)) + 77] defined as the range from the lowest observation to the highest observation with an added extrapolation value 1) to each side. This value can be adjusted depending on the extrapolation capabilities desired from model samples, and how heavy-tailed the data is thought to be. We found 7 = 0.1 to work well, given the data is pre-scaled to have standard deviation 1.
Next, for gridded, equally spaced ordered values { with 3 = min(2 7 and H max( + 7). we evaluate the edf at each of these point. This gives
equally sized set of points {= with ú; == P(x:) E [0, 1}, where K is a hyperparameter guaranteeing the exactness of our approximation with higher value of K in exchange for an increased computational cost. To encompass the complete range of [0, 1}, we set u1 and ux = We fix K 1000 in experiments.
We then use the set of gridded cdf values {ū] 1 to construct an approximation of p-1 through linear interpolation. More specifically, given a value [ū,, uy+1l with 1 S S K 1 that we wish evaluate the inverse to obtain
p-1(u), have
i Z + u -, (ūz+1 - == = f-)
Ux+1 ūj
By consequence, we can easily obtain the gradient of ē w.r.t. marginal hyperparameter p
05 0
u 01
fj Wjt1 Uj
(T - x.))
(13)
0p
0p
(āj+1
-3z)
((ūs++ - -u) Oū, (u ū;)
0ūj+1
(14)
(uj+1
- n;
3p
Op
(Tz+11 - Fzg)
(P(2))1) u) OP(r,)
AP(5f+1)
-(u - po2,))
(15)
(P(2,+1)- - P(x,))²
Op
Op
As such, gradients of the inverse edf become gradients of the edf, which we can efficiently compute with automatic differentiation software.
29

---

Table 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
50
100
300
500
10a
3.
1.73
2.1540.08
0.94
2.43
10
36.47+187
17.14+151
12.82+0.36
7..0+0.22
7.91
30
17.94.1.08
11.04.0.35
10.35
12.87 10.17
9.85.1.40
9T.0917.55
550.5
12.20
48.50
0.73
34.98 10.31
38.59 +4.31
25.82
26.14-0.01
50
115.64 13.06
T12.1642.00
71.43.1.65
78. 78.20
100
268.88
E.1
Comparison to normalising flow on Gaussian Mixture Model
We assess the performance of the QB- -ine mixture of 4 non-isotropic Gaussians across range of dimensions and sample sizes. We simulate n d-dimensional data points from
ply) Σ Th Φ(y: ph: Ex),
where (TI, T2, T3, TT) = (0.2, 0.1, 0.4) and
Ph i.ddd.
ul-50,50]- 50)6,
Σg Lid Wishart Ia).
We compare the Vine with the RQ-NSF as a benchmark off-the-shelf estimator. The hyperparam- eters for the RQ-NSF were chosen to give the best performance on training data, and are 100,000 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 50, and 128 for d = 100. Our results in Table show that the QB-Vine consistently outperforms the RQ-NSF for the dimensions and sample sizes considered.
We additionally considered some even higher dimensional examples to asses the QB-Vine's We study the performance in d 400, 500, 600 dimensions on Gaussian mixture models (GMMs) with 20 random means (drawn uniformly from [-100, 100)d) and non-isotropic covariances drawn from Wishart distribution, with n. 20000 observations, and using 50/50 /50 split for training and testing. We compare the QB-Vine against the RQ-NSF taken as a benchmark for high-dimensional modelling, with the same hyperparameters from the experiments of Table 6 We repeated this study times with different seeds, leading to different GMM models. Figure 2 shows the LPS over the runs for each method and dimension. Further, as the generating distribution is known, we sample from the fitted models and evaluate their samples under the true GMM density (known as the reverse KL divergence, lower better), reported in Figure Finally, we compute the Maximum Mean Discrepancy (MMD)
Vine's 's scalability.
commonly used measure that compares samples to observations, reported Table 7 The MMD assess how close the model is to the true data-generating process by comparing model samples to observed data, with a lower MMD score implying better fit for the data. The results suggest that the QB- Vine has better density estimation as well as sampling capabilities for these examples.
30

---

2991 0002 UMR SAMEE D00
09911 UM613 Lon. 3004 0901
LNNL SNP 3AP
201.
200
0M0 GAR4 1AMS
D0E CA0O 1A0 sANe 2013
09911 1802 MAN DMN
GMM1
CANI
Figure 2: LPS (lower better) for test points on 5 GMMs and for dl 400, 500, 600. The QB-Vine achieves lower LPS values on average than the RQNSF across all 5 GMMs.
Dimension
GMM 1
GMM 2
GMM
GMM4
GMM
29.2213
29.2775
29.0477
29.2993
29.3515
400
29.7847
29.7731
29.7247
29.7835
29.8447
32.7893
32.8354
32.5520
32.6044
32.7249
500
33.1785
33,440T
"33.201
33.4355
33.4443
35.7948
35.8328
35.9390
35.6756
35.7731
600
36.5586
36.6095
36.4700
36.4400
36.7090
Table 7: Comparison of the MMD (lower is better) computed samples from the Vine and RQNSF models across different dimensions and GMMs. Each cell shows the QB Qine value on top and the RQNSF value the bottom, separated by dotted line. The QB- -ine outperforms the RQNSF in cases considered, demonstrating better sample quality via this metric.
GMM 1 dim=400
GMM 1 dim=500
GMM dim=600
10P
Reverse KL under true model
GMM 2 4im=400
GMM 2 dim=500
GMM 2 dim=600
10P
:04
GMM 3 dim=400
GMM 3 dim=500
GMM 3 dim=600
201
101
101
10t
10
GMM dim=400
GMM dim=
GMM dim=600
101
10F
ccn
10h
GMM 5 dim 400
GMM 5 dim=500 500
GMM 5 dim=600
100
10
RQNSF
RONSF mean
QBVine
QBVine mean
Figure 3: Reverse KL divergence (lower is better) for the GMM experiment in high dimensions, assessing sample fidelity. Both models perform equally well.
31

---

NeurIPS Paper Checklist
1. Claims
Question: Do the main claims made in the abstract and introduction accurately reflect the paper's contributions and scope?
Answer: [Yes]
Justification: "We claimed to provide a new density estimation scheme with analytical form which used smaller data to train and can handle large dimensional space. We showed that our proposed QB-Vine has outperformed other benchmark methods for many benchmark dataset the experimental section, providing strength our claim. Further proved theoretically that our proposed method has dimension-indeeendde convergence rate when the dependency of the true model satisfies the assumptions of simplified vine copula.
2. Limitations
Question:
Does the paper discuss the limitations of the work performed by the authors?
Answer: [Yes]
Justification: The main limitation of the method is the of simplified vine copula which is an approximation to the true distribution. We discuss this the discussion but also point that this apparent limitation gives advantages in the sense of fast computation and achieving best results compared to the benchmark methods.
3. Theory Assumptions and Proofs
Question: For each theoretical result, does the paper provide the full set of assumptions and a complete (and correct) proof?
Answer: [Yes]
Justification: We provide the statements of main theorems in Section 3 and the proofs the theorems in the Appendix. To our belief, we have explained all the assumptions (eg. assumption of simplified vine copula) in the statements and we believe the proofs provided are correct and rigorous.
4. Experimental Result Reproducibility
Question: Does the paper fully disclose all the information needed to reproduce the main perimental results the paper to the extent that affects the main claims and/or conclusions of the paper (regardless of whether the code and data are provided not)?
Answer: [Yes]
The experiments are done on publicly available benchmark datasets using other benchmark methods. We provide all the code necessary to reproduce the results as supplementary material.
Justification:
5. Open access to data and code
Question: Does the paper provide open access to the data and code, with sufficient instruc- tions to faithfully reproduce the main experimental results, as described in supplemental material?
Answer: [Yes]
Justification: We have provided the data and code as the supplementary material with this paper. Further host all of this in a personal Github repository, as soon the reviewing is done we would make that repository public for better public access to our code and data.
6. Experimental Setting/Details
Question: Does the paper specify all the training and test details (e.g., data splits, hyper- parameters, how they were chosen, type of optimizer, etc.) necessary to understand the results?
Answer: [Yes]
Justification: In the experimental section we have provided all the details needed for our experiments. In addition, the hyperparameters and their tuning was described in tSection 3. So to our belief, any reader of this paper should be able to implement the algorithms and reproduce exactly the same results following paper.
32

---

7. Experiment Statistical Significance
Question: Does the paper report error bars suitably and correctly defined or other appropriate information about the statistical significance of the experiments?
Answer: [Yes]
Justification: We did provide the LPS as the diagnostics for the performance our method. To show the significance of this method and the compared benchmark ones, have provided means and variances of these key diagnostics over repetitions.
8. Experiments Compute Resources
For each experiment, does the paper provide sufficient information on the com- puter resources (type of compute workers, memory, time of execution) needed to reproduce the experiments?
Question:
Answer: [Yes]
Justification: In the supplementary materials, we provide code files to replicate experiments as well as a notebook file with instructions. In Appendix we describe the compute used and disclose the time taken to run experiments.
9. Code Of Ethics
Question: Does the research conducted in the paper conform, in every respect, with the NeurlPS Code of Ethics
Answer: [Yes]
Justification: This paper provides a new algorithm for density estimation and classifica- tion/regression. As this is a purely an algorithmic research, we believe our work fully satisfies the necessary code of ethics.
10. Broader Impacts
Question: Does the paper discuss both potential positive societal impacts and negative societal impacts of the work performed?
Answer: [Yes]
Justification: This is foundational algorithmic work, which may have potential to impact all sectors of society where density estimation and classifcationn/egressiin is used. We do not believe this proposed methodology has any negative societal impact.
11. Safeguards
Does the paper describe safeguards that have been put in place for responsible release of data models that have high risk for misuse (e.g., pretrained language models, image generators, or scraped datasets)?
Question:
Answer: [NA]
Justification: We only use publicly available data and code which our best knowledge does not constitute high risk for misuse. Hence we do not think this is necessary for our work.
12. Licenses for existing assets
Question: Are the creators or original owners of assets (e.g., code, data, models), used in the paper. properly credited and are the license and terms of use explicitly mentioned and properly respected?
Answer: [Yes]
Justification: All the other data and codes for benchmark methods used in this paper are properly cited and we believe we followed the best practice.
13. New Assets
Question: Are new assets introduced in the paper well documented and is the documentation provided alongside the assets?
Answer: [Yes]
The paper proposes new algorithm, which is thoroughly explained and code is provided.
Justification:
33

---

14. Crowdsourcing and Research with Human Subjects
Question: For crowdsourcing experiments and research with human subjects, does the paper include the full text of instructions given participants and screenshots, if applicable, as well as details about compensation (if any)?
Answer: [NA]
15.
Institutional Review Board (IRB) Approvals or Equivalent for Research with Human Subjects
Does the paper describe potential risks incurred by study participants, whether such risks were disclosed to the subjects, and whether Institutional Review Board (IRB) approvals (or an equivalent approval/review based on the requirements of your country or institution) were obtained?
Question:
Answer: [NA]
Justification: As our studies only involve publicly available datasets not involving any human participation, this approval is necessary for our work.
34


Reviews:
[{"review_id": "0ZVJJDQR5t", "reviewer": "Reviewer_iVqh", "summary": "Inspired by previous quasi-Bayesian (QB) methods, the recursive decomposition of Bayesian predictive posterior distributions, and vine copulas, the work introduces a new adaptation of QB to higher dimensions. The driving idea, introduced in Section 3 consists of making an adaptation of the previous decomposition for univariate densities to a higher number of dimensions, which ends up being a product of marginals times a vine copula with a certain design. Such vine copula models the conditioning among dimensions since the marginals do not. This one works via a decomposition into $d(d-1)/2$ elements that capture the dependency structure. Experiments on UCI datasets with a log-predictive metric show improvement wrt the selected SOTA methods with datasets of small N.", "questions": "**Q1** -- What is the scalability of the simplified vine copula in L146-L152? Could the authors add some extra info about the computational cost per iteration n?\n\n**Q2** -- Why so much focus on the log-predictive metric? Wasn't there an additional informative metric that could give the reader more insights about the performance compared with other inference methods? Is there a log() operator missing in the equation on L271?\n\n**After rebuttal comment:** Dear authors, thanks for your rebuttal and the clarifications made regarding my points of weakness and concerns. Some things are clearer now to me, so I am happy to increase my score to 5 now. However, I still feel that some issues raised in my review have not or at difficult to fix atm, and I feel that there are not yet strong reasons to clearly accept the paper (i.e. clarity/readability and empirical results are yet limited, and I still have doubts on the novelty side). I don't object to acceptance, but won't be the supporter of that either. Best of luck.", "limitations": "Yes", "rating": 5, "confidence": 3, "soundness": 2, "presentation": 3, "contribution": 2, "strengths": "I would like to add some points of strength that I think are worth to be mentioned:\n\n- The work introduces well enough the problem of higher dimension d in certain Bayesian inference problems and how some SOTA methods and frameworks might struggle when this is larger than a few dozens. Additionally, the reference to previous works and advances is important and well pointed out in my opinion, giving the right credit for the technical parts to each one of the previously published works.\n\n- The paper is in general concise, and I could follow the technical details -- so I don't think there are details missing in the math part, despite some lack of analysis in certain directions that I'd had loved to see.\n\n- In general, the idea of breaking the conditioning in this way, exploiting both recursive predictives and an additional element that captures correlations in a moderately-scalable way (i.e. vine copulas) is interesting to me and for the paper.", "weaknesses": "Some comments that I would like to add in terms of weak parts or ideas that I consider kind of a problem (or at least I'm concerned about):\n\n- I think the technical derivations, and in general the notation could be improved. From the usual Bayesian perspective and for readers familiar with probabilistic methods, it is not really orthodox and the subscript-superscript system could be confusing at times. The work introduces too many things sometimes instead of focusing on clarity and highlighting the actual contribution and strong ideas proposed in the manuscript. At least, that's what I get from reading it.\n\n- The contribution and novelty are a big concern to me. Honestly, I feel that if Eq. (7) is the big contribution or the point of novelty of the paper, it would not be enough, as it is just an adaptation of the previous derivations in Section 2 / Background to a vine copula. Additionally, the ratio in Eq. (8) is trivial, right? nothing really to derive there, or am I wrong here? \n\n- The experiments are very limited, on small UCI datasets and only show one empirical result which is the improvement with the dimensions, however, no discussion about complexity/computational cost, run times, sensitivity to the choice of the vine copula system, etc is made. In that regard, I fear a bit how scalable really is the vine copula with such $d/(d-1)/2$ mechanism to capture conditioning.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "# Reply to weakness 1\nWe appreciate your recommendation on this. Following your and other reviewer's comments, we have (1) interchanged subscripts and superscripts between dimension and predictive step, (2) added an introduction of vines as a subsection in Section 2, providing a clearer overview there and subsequently focusing more on our contribution on Section 3, along with minor fixes.\n\n# Reply to weakness 2\nThe big contribution of our paper is the construction of a recursive QB method that is not restricted by assumptions on the DPMM kernel. As explained in point 1 of our main reply, the main drawback of previous methods is their necessary assumptions on the multivariate structure to obtain their multivariate recursion, leading to sub-optimal fits on multivariate data. Our tool to avoid restrictive assumptions on the multivariate structure is Sklar's theorem, which in itself is not new, but its application to bypass kernel assumption in a DPMM model is novel. Previous work [31] made incremental improvements to the recursive construction, whereas here we take a wholly different approach, overcoming the limitations of the existing literature in the field, and outperforming them. Additional points of novelty are the construction of a recursive estimator that can be computed in parallel instead of sequentially like the auto-regressive decomposition of [25,31], allowing multiple-fold computational savings, which is the primary focus of this line of research for the DPMM. This is further assisted  by optimising the energy score, halving our training time compared to existing work. Equation (8) is included there to make the potential for parallelisation and unrestricted dependence structure apparent.\n\nOur approach can further benefit from improvements to marginal recursion (as e.g. in [29] or future works since this is an active area of research as pointed out by reviewer qqFA) as well as improvement to copula models since both the R-BP and vines can be replaced with upgraded methods. This places our work in a unique intersection of recursive computations and multivariate copula models, connecting these up to now mostly disjoint subfields and opening a flow of cross-fertilisation between them. \n\nWe thank you for your comment and take this as our queue to better motivate our contribution in the introduction. For the camera ready, we will add more context about the DPMM assumption (following point 1 in the main reply) in the text to clarify our innovation.\n\n\n# Reply to weakness 3:\nWe appreciate your concern about scalability. We have included an example on Gaussian mixture models in high dimensions addressing this (see main reply). Further, to asses sensitivity, we expanded our study on the digits dataset, totalling 15 runs per each training size. We report that the change of copula structure does not negatively affect the final predictions, demonstrating robustness to the vine structure. We included a pseudo-code of our proposed algorithm in the pdf page.\n\nThe computational cost of each step mirrors that of the existing methods, R-BP and vines [65], with the exception that we half the time of the R-BP [25,39] due to the optimisation of the Energy score. For the R-BP, our decomposition into marginals yields the same cost as the cheapest initial update step of [25] or [31] across all our recursive steps. Through parallelising, we obtain a constant cost with $d$, and scale as in [25,39], meaning $\\mathcal{O}(n^2)$ for initialising the recursion (the first bullet point of the algorithm in the pdf), and $\\mathcal{O}(n)$ to compute the pdf or cdf at a point. For the vine, the number of pair copulas grows quadratically with the dimension, and the number of tree structures is exponential with the dimension, leading to greedy algorithms being used for the tree selection, see [66]. To simulate from a vine, one can efficiently compute the inverse probability transformations (see e.g. Chapter 6 in [16]). Further, one can \"truncate\" the tree by assuming that pair copulas past a certain degree of conditioning are uniform, reducing the complexity significantly. One can also set a threshold of a simple association measure like Kendal's Tau, computable without model fitting, to decide which pair copulas to model and which to keep a priori. The drawbacks of the vine can be mitigated with truncation or thresholding as discussed above but inevitably loses out on modelling power. The flexibility of the vine in that regard has been studied in the literature, we refer to [16,17,65,66,67,87] among others. Vines are generally accepted as the best copula model for high dimensions, hence our use of it. \n\n\n# Reply Q1: \nThank you for this important point. Following your question, we included an algorithmic description of our method (see reply to all reviewers) as well as a discussion on computational times and complexity.\n\n# Reply Q2:\nThe LPS is a strictly proper scoring rule, meaning that the unique minimiser of the LPS is the true data generating mechanism. Further, the LPS is equivalent to the negative log-likelihood as commonly used to evaluate density estimators in the literature. Notably, the LPS is the only metric used in the previous QB literature [25,31]. As we compare with the results of [31], and their code has been removed from their website, we are unfortunately unable to replicate their experiments and evaluate other metrics. Following your insight, in our additional high-dimensional experiment on Gaussian Mixture models, we also compare the sampling quality of the model by computing the density of the samples under the model density (known as the reverse KL divergence) as well as the MMD, observing similar performance. There was indeed a missing log() on L271, we thank you for spotting it and have corrected it.", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}, {"review_id": "q2yISACEvJ", "reviewer": "Reviewer_QJH9", "summary": "The authors propose a novel method for modeling high-dimensional distributions (for density estimation and supervised learning), where they break the estimation task into estimation of univariate marginals and estimation of a multivariate copula. To expedite the univariate estimation tasks they utilize the novel quasi-Bayesian (QB) estimation method, and they use a simplified vine copula to approximate the multivariate copula. They present empirical results that compare their method to others in density estimation and supervised learning.", "questions": "- L29: The paper makes a confusing start regarding its notation. $n$ is used here without being introduced, on L69 it's referred to as $K$, which leads to further confusion. This is likely the most important variable in the paper, so a proper introduction is required - especially because the recursive nature of the approach in question may not be as obvious for a reader from a different subfield.\n- L163: Potential implications of using simplified vine copula can be discussed here, with a reference to forthcoming Thm 3.3.\n- L174: Any particular reasons for the choice of Cauchy as the initial distribution?\n- L188: Please use the additional space afforded after reviews to make this discussion more explicit.", "limitations": "The authors provide a satisfactory discussion of the limitations of the paper.", "rating": 7, "confidence": 2, "soundness": 3, "presentation": 3, "contribution": 3, "strengths": "- The authors present a well-defined problem that is of potential interest to conference readership, place it within the existing literature, and clearly demarcate their innovations.\n- They break down the estimation problem into two sub tasks (estimation of the univariate marginals and the copula), and offer tangible contributions for both tasks, and provide convergence results for their estimator.\n- Empirical results are convincing regarding authors' method's practical usefulness.", "weaknesses": "Although mostly outside of my expertise, this is a well-written paper overall. However, there are a number of potential changes that can improve the readability and accessibility of the paper. I list these in the section below.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "We appreciate your effort in reviewing our paper and are glad you found it well-written - thank you. We provide answers to your suggestions below. \n\n# Reply Q1:\n We thank the reviewer for pointing out this important oversight and have accordingly removed all $K$ and replaced them with $(n)$. We also removed the notation in the introduction and changed the indexing to be more standard, with dimensions as superscripts and recursive numbering as a subscript. \n\n# Reply Q2:\nWe thank the reviewer for their comment. Following other reviewers' suggestions, for clarity's sake, we have moved the introduction of vine copulas to Section 2 as a subsection. We will incorporate your remarks there, detailing the implications of a simplified vine copula.\n\n# Reply Q3:\nThe choice of initial predictive has 2 implications. Firstly, it is an implicit statement about your beliefs on observables, in a similar way to a prior. Secondly, it contributes to the efficiency of the recursion in fitting data. We discuss both of these aspects in Appendix E. In particular, a Cauchy distribution is effective at minimising numerical overflow when computing the cdf of data in the tails, making our algorithm more robust to different types of data. Further, taking a heavy tailed distribution such as the Cauchy coincides with the theoretical work on similar recursive density estimators, where it is common to assume that the tails of the true data generating density are not too heavy compared to the tails of the predictive, see e.g. condition (3.2) in [59].\n\n# Reply Q4:\nWe thank you for this comment. We will expand upon this point in the camera ready.", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}, {"review_id": "UBU5dE6XTi", "reviewer": "Reviewer_2HTp", "summary": "This paper develops the recently proposed quasi-Bayesian methods by applying vine copula (hence named QB-Vine) to the recursive Bayesian predicative distributions and bypassing the need for expensive posterior integration. The proposed method consists of two parts:  independent recursion of marginals by bivariate Gaussian copulas; and estimating the simplified vine copula to capture data dependence hence relaxing the kernel assumption for DPMM. Error bounds for both the distribution functions and the copulas are provided to justify the proposed QB-Vine. Numerical examples including density estimation and regression/classification are used to showcase the advantage compared with the state-of-the-art alternatives.", "questions": "Should there be $|\\mathcal S_{ij}$ after the second conditional distribution in equation (10)?\nDid you have dimension specific bandwidth $b$ and correlation $\\rho$ for all the numerical experiments? Did it increase the overall computation time significantly by having them different for each dimension?\nLine 223: \"convergence\" should be \"converges\"\nFigure 2: why LPS increases from 500 samples to 900 samples for QB-Vine?", "limitations": "See weakness", "rating": 6, "confidence": 4, "soundness": 3, "presentation": 3, "contribution": 2, "strengths": "The paper proposes a novel Bayesian method to compute the predictive distribution. The proposed method, QB-Vine demonstrates numerical advantage over alternatives. Theoretic characterization on the errors is provided.", "weaknesses": "The dimension of the problems is relatively low (up to 64).", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "We thank the reviewer for their effort in reviewing and offer answers to their comments below.\n\n# Reply to weaknesses: \nTo showcase the scalability of our approach in higher dimensions, we have expanded the experiments from Table 6 in the paper with a study in $d=400,500,600$ dimensions on Gaussian mixture models (GMMs) with 20 random means (drawn uniformly from $[-100,100]^d$) and non-isotropic covariances drawn from a Wishart distribution, with $n=20000$ observations, and using a $50/50$ split for training and testing. We compare the QB-Vine against the RQ-NSF taken to be a benchmark for high-dimensional modelling, with the same hyperparameters from the experiments of Table 6. We repeated this study 5 times with different seeds, leading to different GMM models. In the additional uploaded material, we included Figure 3 showing the LPS over the 5 runs for each method and dimension. Further, as the generating distribution is know, we sample from the fitted models and evaluate their samples under the true GMM density (known as the reverse KL divergence, lower is better), reported in Figure 1. Finally, we compute the MMD between samples and observations, reported in the table below.\n\nDimension  GMM 1      GMM 2      GMM 3   GMM 4   GMM 5\n-------------------------------------------------------\n--------------------------------------------------------\n400\n------------------------------------------------------\n--------------------------------------------------------\nQB-Vine 29.2213    29.7847    29.2775 29.7731 29.0477 \n------------------------------------------------------\n--------------------------------------------------------\nRQNSF 29.7247    29.2993    29.7835 29.3515 29.8447 \n-------------------------------------------------------\n--------------------------------------------------------\n500\n------------------------------------------------------\n--------------------------------------------------------\nQB-Vine        32.7893    33.1789    32.8354 33.4401 32.5520\n-------------------------------------------------------\n--------------------------------------------------------\nRQNSF           33.2011    32.6044    33.4355 32.7249 33.4143 \n-------------------------------------------------------\n--------------------------------------------------------\n600\n------------------------------------------------------\n--------------------------------------------------------\nQB-Vine        35.7948    36.5586    35.8328 36.6095 35.9390\n-------------------------------------------------------\n--------------------------------------------------------\nRQNSF        36.4700    35.6756    36.4400 35.7731 36.7090 \n-------------------------------------------------------\n--------------------------------------------------------\n\\caption{Comparison of the MMD (lower is better) computed on samples from the QBVine and RQNSF models across different dimensions and GMMs. Each cell shows the QBVine value on top and the RQNSF value on the bottom, separated by a dotted line.}\n\n# Reply Q1: \nYes, this is a consequence of the vine decomposition following Equation (9), meaning the copula density $c_{i,j}$ is itself conditional on $S_{ij}$. The vine copula of equation (10) is still a full vine copula, with no simplifying assumption. The simplifying assumption of a simplified vine copula then ignores that last conditioning in favour of a more parsimonious model. Additionally, the simplifying assumption makes the vine model easier to estimate, as then all the pairs $(P_{i|S_{ij}}(x_{i|S_{ij}}),P_{j|S_{ij}}(x_{j|S_{ij}}))$ can be used in the estimation of $c_{ij}$ instead of having to also take into account the values of the conditioning set $S_{ij}$. We take your comment as a potential confusion induced by our sentence on L154 \"rewriting any conditional densities as copulas\", and have therefore removed it for clarity.\n\n# Reply Q2 \nWe use a dimension-dependent correlation $\\rho$ (L176), but use the same bandwidth for all KDE pair copulas (L181). We have modified L176 to state this more clearly.\n\n # Reply Q3\n The computational time is the same as estimating a common correlation parameter across all dimensions, since each recursion still has to be computed to select $\\rho$. In fact, choosing a separate $\\rho_i$ per dimension is preferred since it allows to estimate the parameters in parallel without each dimension needing to interact, for example to pool gradients as would be the case with a common bandwidth. Therefore, with parallelisation, selecting a different bandwidth for each dimension takes the same time as selecting the correlation for a single dimension.\n\n# Reply Q3 \nWe thank the reviewer for spotting this typo and have fixed it.\n\n# Reply Q4 \nWe believe this was a consequence of not taking enough runs in the experiment. We have now updated our picture to include 15 runs, showing a smooth decrease of the LPS with training size, as one would expect.", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}, {"review_id": "qTqs8ESuWP", "reviewer": "Reviewer_qqFA", "summary": "Quasi-Bayes, in the sense of building a model as a consistent sequence of predictive distributions, has recently received a lot of interest. The initial construction in [25] used products of bivariate copulas to build their predictives; the authors propose here to extend the construction to vine copulas, a flexible model for multivariate distributions.", "questions": "Q1. Can you give pseudocode for your method applied to a particular task, including all hyperparameter tuning steps like pair-copula estimation? \n\nQ2. Can you then discuss the computational cost of each step, and identify the bottlenecks?\n\nQ3. A convincing answer to items 11-12-13 in my list of major comments would also make me reconsider my mark.", "limitations": "This is fundamental work, and no immediate negative societal consequence is foreseen.", "rating": 5, "confidence": 4, "soundness": 3, "presentation": 2, "contribution": 3, "strengths": "* Quasi-Bayes in the wake of [25] is indeed a hot topic in Bayesian statistics, and the paper is potentially impactful.\n* The initial copula constructions in [25] felt a bit constrained, and more flexibility would intuitively help.\n* I overall enjoyed reading the paper, though I have many (hopefully constructive) comments below.", "weaknesses": "## Major\n1. p1 L1 I would avoid flattering formulations like \"[25] initiated a new era of Bayesian computation\" and L33 \"[25] heralded a revolution in Bayesian methods\", \"liberating Bayes\". While you are free to think it and predict a revolution, that kind of statements about a 2023 paper are necessarily statements of opinion or personal predictions at best, and we should let time tell whether a new era or a revolution happened. You can say, e.g., that [25] proposed a stimulating change of paradigm, that is a fact. More pragmatically, starting a paper with opinion statements puts the reader in a suspicious mood.\n2. Abstract: the sentence \"But extensions [...] used\" is unclear. What assumptions? What kernel?\n3. Figure 1 is not particularly useful, I'd say, and the fonts in the decomposition are a bit strange. I know ICLR papers, for example, tend to have a small summary figure at the beginning, but I believe that it doesn't especially help here. The space gained could be used to give e.g. a more formal introduction to vine copulas in the main text.\n4. p2 L42 \"DPMM kernel structure\". Maybe remind the reader on p1, when you first mention the DPMM, of what it is and what you call the kernel. It will not be obvious to many.\n5. The introduction should contain a sentence defining vine copulas, or \"vines\". It is likely not a standard notion in the NeurIPS community, and the introduction should make it clear what the title of the paper means. Actually, the definition is somehow informally given in the caption of Figure 1, but I had missed it at first reading. I would avoid putting anything else in a caption than the description of what the figure shows. Important details should go to the main text.\n6. p2 L56 \"under the assumption of well-identified simplified vine copula model\". It is unclear at this stage what you mean here. \n7. Eqn 2: as far as I understand, going from the first to the second line requires Sklar's formula, at least if I follow the derivation of [25]. If you also use it, this should be mentioned. There is a reference to Appendix A, where copula tools are introduced, such as Sklar's theorem, but I'm not sure the argument is explained.\n8. p3 L 82-82, in what sense does $c^{(n)}(\\cdot,\\cdot)$ converge to $1$? How does is guarantee a convergence for $p^{(n)}$? Why is the copula symmetric?\n9. The paragraph L102--L113 is a bit dense and can gain in clarity. In what sense and under what assumptions does the \"univariate R-BP model converge to a limiting distribution ... with density...\"? What do you mean by \"Newton's algorithm\"? Also, there are steps to explain if you want to talk about how quasi-Bayes \"approximates a posterior density\"; so far, we have only talked about prediction. What parameter do you consider? How do you approximate the posterior density? You could refer to the corresponding passages in [25] and explain them in a few words. Or avoid discussing parameter inference if it is not central to your contribution.\n10. p4 L128 what do you mean by \"the predictives are unconditional marginal densities\"? The sentence is puzzling.\n11. Importantly, the authors of [25] restrict their designs of predictive distributions so that they satisfy their martingale condition (Eqn (4.1) in the Arxiv v2 version). They use their Corollary 1 for univariate copulas, and note in Section 4.3 that a multivariate extension of their Corollary 1 is not easy, which prompts them to take instpiration from  factorized kernels in DPMMs, for which they can guarantee (4.1). So the obvious question is: how do you guarantee that your vine predictives satisfy Fong el al.'s martingale condition? Or, alternatively, how do you guarantee the existence of the limiting martingale posterior $P^{(\\infty)}$? It is possible that I am missing an obvious argument. \n12. Lemma 3.2 uses $P^{(\\infty)}$, which should be formally defined. Or at least, a reference should be given to the formal definition in [25]. But I think that giving a formal definition would help, since it would entail checking the martingale assumption of [25]; see bullet 11 above.\n13. Overall, guarantees like Lemma 3.2 or Theorem 3.3 are guarantees in probability, under the distribution formed by the product of your predictives, am I right? In that case, these guarantees are meaningful only if I believe that the product of predictives is actually a good model for the data generating process. Is this realistic? Even under the classical simplifying assumption on conditional copulas that you are making?\n14. p5 L146-173 the introduction to vine copulas is a bit confusing and informal. This section should be more formal, as understanding vine copulas is key to the paper, and I think not standard to many NeurIPS readers. I suggest having a subsection of Section 2 that formally introduces vine copulas, independently of their application to quasi-Bayes, and gives your example (23). Then Section 3 can focus on how you use vine copulas to define your joint predictive distributions. Also, I would use Pseudocode at the end of Section 3, to summarize all the steps in your procedure: hyperparameter tuning, estimation of the pair copulas, etc. I am not 100% confident that I can write this pseudocode from the text only.\n15. Relatedly, having pseudocode would allow you commenting on the computational cost of the different steps. p9 L335 you mention \"a search over a vast model space during estimation\", but this has not been commented on before. \n16. Theorem 3.3: what do you mean by \"Assuming a correctly identified simplified vine strucgture for $c^{(\\infty)}(\\mathbf{u})$?\". If I understand well, a more standard phrasing would be to formally introduce your model, say as $P$, and then state that under the distribution $P$, such and such is true. \n17. The caption of Table 1 mentions error bars of two standard deviations averaged over 5 runs. I am unconfortable with the number of runs, which is intuitively too low for a $2\\\\hat{\\sigma}$ interval to make sense, be it through a CLT or a Chebyshev argument.\n\n## Minor\n* p1 L29 the notation $\\mathbf{p}^{n}(\\mathbf{x})$ has not been defined yet, so I would simply avoid using the notation here.\n* p2 L52 while you are not the only ones to do so, I don't see the need to boldface \"main contributions\".\n* Theorem A.1: I would define what a copula distribution is.\n* p3 L87 \"intractable\": following [25], I would even say it is not particularly desireable, as they try to bypass the need to specify a prior-likelihood pair. \n* In Eqn (8), the densities should appear with their arguments. \n* Section 5: don't you miss a log in your expressions of the log score?", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "# Reply 1:\nWe did not mean to alienate readers, and will modify the text according to your suggestions. For the camera ready, we will use your phrasing for the abstract and use more neutral expressions in the introduction, citing papers emanating from [25] as evidence of an active area of research instead. \n\nReply 2:\nPlease see point 1 in the main reply.  \n# Reply 3,5,6,14:\nWe appreciate your judgement and removed the figure. We use the space to describe vines briefly in the introduction and thoroughly as a subsection of Section 2, accordingly correcting the relevant passages in the rest of the text.\n\n# Reply 4:\n We thank the reviewer for this suggestion and have implemented it in the text in the style of point 1 of the main reply.\n\n# Reply 7:\n We have added a sentence below the equation explaining the use of Sklar's theorem in the derivation by the definition of a copula.\n\n# Reply 8:\nIn Corollary 1 of [25], for univariate predictives, it is shown that the existence of a copula sequence for predictive recursion is equivalent to the sequence of predictive densities satisfying the martingale condition, in turn implying that $p^{(n)}(x)\\rightarrow p^{(\\infty)}$ almost surely for each $x\\in\\mathbb{R}$ by the martingale convergence theorem. Due to the recursion of Equation (2), this implies $\\lim_{n\\rightarrow\\infty}c^{(n)}(x,x^{n})=1\\, a.s. \\forall x$, meaning that the copula does not interfere with the convergence, and eventually the predictive does not change anymore. We have made the almost sure convergence explicit in writing and have justified it by saying \"...$\\lim_{n\\rightarrow\\infty}c^{(n)}(x,x^{n})=1\\, a.s. \\forall x$ as a consequence of the almost sure convergence of $ p^{(n)}$ with $n$ [25]\". Lastly, the copula is symmetric as one can freely exchange the components of the copula in equation (2). Namely, exchange $f(x|\\theta)$ with $f(x^n|\\theta)$ for the expression of the joint density in the numerator and $p^{(n-1)}(x^n|x^{1:n-1})$ with $p^{(n-1)}(x|x^{1:n-1})$ for the marginals in the denominator, thus obtaining $c(x_1,x_2):=\\frac{p(x_1,x_2)}{p_1(x_1)p_2(x_2)}=\\frac{p(x_2,x_1)}{p_2(x_2)p_1(x_1)}=:c(x_2,x_1)$.\n\n# Reply 9:\n We are quoting the result from theorem 5 in [25], saying that the R-BP distribution $P^{(n)}$ converges in total variation to $P^{\\infty}$ which is absolutely continuous with respect to the Lebesgue measure with an associated density $p^{\\infty}$. The theorem needs the R-BP density $p^{(n)}$ to be continuous, $\\int_K {p^{(n)}}^2(x) d x<\\infty$ for $K$ a compact subset of $\\mathbb{R}$ with finite Lebesgue measure, the weight sequence to have the form\n$\\alpha_i=\\left(2-\\frac{1}{i}\\right) \\frac{1}{i+1}$ and that the parameter of the copula be $\\rho<1 / \\sqrt{3}$. Then, in [39] the R-BP density is shown to converge to the true density in a Kulbach-Leibler sense, under similar assumptions. Newton's algorithm is the original work that initiated the study of recursive density estimates for the DPMM. The two cited papers [68,69] focus on establishing a recursion for the mixing distribution of a DPMM. However, at every step of the recursion their approach needs to solve an integral with dimensionality scaling with the dimension of the data, thereby making it impractical in higher dimensions. A thorough review of the original and ensuing works has been provided in [59]. We comment on this in Section 4. The \"posterior density\" was a typo, as we meant the predictive density, which is the object of the R-BP recursion. Parameter inference is indeed not central to our work, so we do not discuss it. We thank the reviewer for commenting on this and have ameliorated the clarity of the paragraph by referring to the corresponding results explicitly and fixing the typo about the posterior.\n\n# Reply 10:\nThe goal was to emphasise that the univariate densities are not conditional as in [25,39], but are independent of each other allowing them to be modelled separately (i.e.~in parallel). We appreciate this comment on clarity and have rephrased the sentence into \"the predictives are simple univariate densities\".\n\n# Reply 11:\nIn the interest of space, we include this as a comment below.\n\n# Reply 13: \nYou are correct, our proofs are statements about the distribution formed by our model of the joint predictive, being marginal R-BPs multiplied with the vine copula. We do believe this is indeed a useful model with a rich complexity. Firstly, the marginal R-BP minimises the KL to the true data generating process with $n$ [39], and converges in total variation [25], leading to precise pseudo-observations $\\{u_i=P^{(n)}(x_i)\\}_{i=1}^d$ (lemma 3.2). Then, we employ copulas to recover the joint predictive, knowing that such a copula exists and is unique (theorem 3.1). We believe the best copula model currently available in the literature is a simplified vine copula, as it is both computationally feasible and offers a large model flexibility with non-parametric pair copulas (see e.g. \\cite{czado2019analyzing}). Further, the convergence rates of the R-BP (lemma 3.2) and the vine copula (theorem 3.3) make them appealing density estimators.\n\n# Reply 14:\n Please see the pseudo code in the pdf, and our point 2 in the main reply.\n\n# Reply 15:\n We will include a discussion on the computational cost of our approach. \n-see reply to other reviewer and experiment on computing time-\\\\\nThe model space we refer to is the specific decomposition of the vine copula, or equivalently its tree structure. We discuss this in L185-187. Following your comment, we will rephrase the sentence in Section 6 to be coherent with our discussion on computation costs mentioned above.\n\nFurther replies and the proof are in the comment.", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}]

Author rebuttal:
#### We thank all the reviewers for dedicating their time and effort to the review of our paper. Below, we outline the main points raised by reviewers and how we address them. We also uploaded a pdf file with extra results and a figure showing the algorithm for the QB-Vine. Further, we have posted an individual reply to each reviewer discussing their specific comments.

# Introduction to DPMM and main novelty of our approach:

Reviewers iVqh and qqFA have asked about the contribution of our work, and clarifications about statements related to previous work. We provide details here and have reformulated relevant parts of the abstract and introduction to clarify the novelty further.
Our paper addresses the recursive modelling of the DPMM's predictive density, also studied by previous work [25,31]. The multivariate DPMM is formulated as
$$f(\mathbf{x} \mid G)=\int_{\Theta} K(\mathbf{x} \mid \boldsymbol{\theta}) d G(\boldsymbol{\theta}), \text { with } G \sim \operatorname{DP}\left(c, G_0\right)$$
where $K$ is a kernel for the observables $\mathbf{x}\in\mathbb{R}^d$ parameterised by $\theta$, similarly to the kernel in kernel density estimation, and $G$ is called the mixing distribution, upon which a Dirichlet process prior is placed with base measure $G_0$ and concentration parameter $c>0$. The specification of $K$ and $G_0$ are thus important as they regulate the expressivity of the model, but their computation in a Bayesian model requires significant effort. To address this shortcoming, [25] provide a copula-based recursion obtained by assuming $K(\mathbf{x}|\theta)=\prod_{j=1}^d\mathcal{N}(x^j|\theta^j,1)$ and $G_0(\boldsymbol{\theta})=\prod_{j=1}^d\mathcal{N}(\theta^j|0,\tau^{-1})$, $\tau>0$, meaning both the kernel and base measure are assumed independent across dimensions, lacking the expressivity required to capture dependencies in the data. In [31], the form of the kernel is relaxed to be autoregressive with $K(\mathbf{x}|\theta)=\prod_{j=1}^d\mathcal{N}(x^j|\theta^j(\mathbf{x}^{1:j-1}),1)$ where the kernel mean $\theta^j:\mathbb{R}^{j-1}\mapsto\mathbb{R}$ is dependent on previous dimensions, and the base measure of [31] is a product of Gaussian Process priors $G_0(\boldsymbol{\theta})=\prod_{j=1}^d GP({\theta}^j|0,\tau^{-1}k)$ for $k:\mathbb{R}^{j-1}\times\mathbb{R}^{j-1}\mapsto\mathbb{R}$ a covariance function.

In our work, we posit the existing forms of $K$ and $G_0$, necessary for their recursions, lack the flexibility to model multivariate data accurately. As a result, we introduce the QB-Vine as a more general approach circumventing assumptions on the recursive form. Instead of deriving an autoregressive recursion through assumptions on $K$ and $G_0$, we show that by applying Sklar's theorem on the joint density, one only needs to specify the recursive form of the marginals, overcoming the limitations of existing work and leading to better performance in experiments. Lastly, our proposed method is inherently parallelizable over dimensions instead of sequential, leading to significant computational gains.

# Scalability to high dimensional problems:

 To showcase the scalability of our approach in higher dimensions, we have expanded the experiments from Table 6 in the paper with a study in $d=400,500,600$ dimensions on Gaussian mixture models (GMMs). We compare the QB-Vine against the RQ-NSF as a benchmark, repeating the experiment 5 times, and showing our competitive performance. In the uploaded pdf page, we show: the LPS (Figure 3), the reverse KL under the true model (Figure 1) and the MMD on samples (Table 1, see the reply to reviewer 2HTp). Further details are in the reply to weakness 1 of reviewer 2HTp.

# Computational Complexity:  
We added pseudo code to the pdf for the QB-Vine. $M$ is the number of permutations, $B_1,...,B_l$ the bandwidth grid, $V$ the cross validation percentage and $J$ the number of samples to compute the Energy score in the vine. Please see our reply to Q1 of reviewer qqFA for more details. We further ran 10 additional runs on the digits dataset in Figure 4 of the pdf, to assess sensitivity. We discuss computational costs in detail in the reply to weakness 3 of reviewer iVqh. 


# Checking the Martingale condition:
Reviewer qqFA asked to prove the martingale condition (4.1) in [25] for the QB-Vine, please see the proof in our reply 11 to them. The only assumption is the ratio of two consecutive vines must be 1, i.e., the dependence structure is constant between predictive steps. We note that these derivations hold for any marginal recursive construction of the form (2) and any copula density used for $\mathbf{c}^{(i)}$, but we focus here on the QB-Vine. We interpret the condition of having the same dependence structure between steps as natural when data is supposed to come from the same data-generating process, which is indeed the circumstance in which we apply the QB-Vine. Further, given $n$ observations, the best guess of the multivariate copula is given by fitting it at the last iteration, which is our approach in experiments.

Meta-review: The authors tackle the problem of high-dimensional density estimation from modest amounts of data via a combination of two techniques: quasi-Bayesian inference, in which Bayesian updates are viewed in terms of copulas (pairing past data and a new observation) that are then approximated in a recursive process; and vine copulas (specifically, simplified vine copulas) in which multivariate distributions are assembled from bivariate copulas (Gaussian, in this paper).  The authors demonstrate superior performance compared to alternative methods on several problems with high dimension (in this case meaning 8-64 dimensional real data, 100-dimensional for one synthetic example in the appendices) and modest data relative to the dimension.  In the rebuttal period, the authors provided a comparison for an even higher dimension synthetic example (up to 600 dimensions).

Though surely more efficient than Monte Carlo-based alternatives, the inference procedure still involves an expensive search over a large model space (the space of possible nested trees is very large, and there are other hyper-parameters as well).  I note that while questions about scalability with dimension were surfaced during the discussion, there were few questions about scaling with data size -- though the authors point out that the one UCI test case where their method was beaten by a competitor involved a relatively large amount of data relative to the dimension.

The reviewers on this paper were weakly positive, though several expressed low confidence in their evaluation.  The most positive reviewers were the least confidence.  The other three landed on two weak accepts and one borderline.  Reviewers liked the combination of ideas, but were nervous about the utility for higher dimensions and data sizes, even with the updates in the rebuttals.  There was also some uncertainty about the novelty in the initial reviews, though this seems to have been largely resolved in the discussion period.

In general, machine learning audiences tend to favor methods that scale well to high dimensions and large data sets.  This paper focuses much more on what is large (in dimension and data size) with respect to the statistics literature, and this seems to be reflected in the reviewer reactions.  However, this is still an important regime, and this is an interesting combination of ideas, and the reviewers (and the AC, in my own reading) seem to agree with this.  For this reason, I am recommending acceptance.
Decision: Accept (poster)