Reply 12:

We note that $P^\infty(x)$ used in lemma 3.2 is univariate, with the martingale condition being established for the univariate R-BP in [25] already. Following your advice, in the camera ready, we will dedicate a paragraph to the introduction of the martingale condition and the proof of the QB-Vine satisfying it.

Reply 16:

We thank you for this comment, and have rewritten the statement with that structure.

Reply 17:

Please see our reply to Q2 of reviewer iVqh, explaining our limitation due to [31]. With their code unavailable, we assume the best case scenario for other methods, believing those intervals are adequate there, but study the intervals of the QB-Vine in more detail, running 15 more runs for each point Figure 2 (see pdf).

Reply minor points:

We incorporated the changes in our manuscript.

Reply Q1:

Please find the pseudo-code in Figure 2 of the pdf, including hyperparameter estimation. Some brief clarification: the training data can be permuted in full or partially depending on how much data is available, amd is captured by the variable $M$ giving the number of permutations. We keep $M=10$. Other hyperparameters are the number of points $B_1,\ldots,B_l$ to select the bandwidth of the vine. Our experiments were run with a grid of size $50$ between $2-4$ for UCI dataset and $0.5-3$ for the rest. $V=0.8$ to have $5-$fold cross-validation, and $J=100$. $\rho$ can be optimised with gradients and converged very quickly in our experiments, requiring less than five evaluations on average.

Reply Q2:

The computational cost of each step comes from the existing method, R-BP or vines [65,], with the exception that we half the time of the R-BP [25,39] due to the optimisation of the Energy score. Due to space constraints, please see our reply to weakness 3 of reviewer iVqh.

Finally, we thank the reviewer many times for their great suggestions and critical comments. We truly believe the paper is stronger as a result, and express our greatest thanks for dedicating your time to our work. We hope we addressed your concerns, and remain available to discuss any further points you would wish to raise.

Reply 11:

We do not claim to converge to a multivariate $P^\infty(\mathbf{x})$. The $P^\infty(\mathbf{x})$ is only relevant when one wants to do predictive resampling to impute the unobserved data, as done in the martingale posterior for parameter inference. We do not investigate predictive resampling and do not advertise our work as such, nor focus on parameter inference. We need the marginals to converge (lemma 3.2). Then, for our theorem 3.3, this is purely a statement on the density as done in the vine copula literature, not on the limiting $P^\infty(\mathbf{x})$ martingale posterior in multiple dimensions. Similarly, theorem 3.1 adapts Sklar's Theorem to multivariate densities, holding for any $n$. For experiments, we show that our model approximates well all the datasets, doing density estimation, for which a $P^\infty(\mathbf{x})$ is not required. Our paper focuses on modelling predictive densities in an efficient way, which is a worthwhile pursuit as classical density estimators such as KDE are known to scale poorly with dimension. However, we were able to prove the martingale condition, see point 4 in the main reply, and the proof below. (We could not get it to format properly, we express our deepest apologies and hope it is still readable.) We write $\mathbf{c}^{(i)}= \mathbf{v}^{(i)}$ for clarity and denote inputs to copulas as vectors $(P_1(x_1),\ldots,P_d(x_d))=[P_i(x_i)]$:

Reply 1:

We 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.

Reply 2: Please see point 1 in the main reply.

Reply 3,5,6,14:

We 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.

Reply 4:

We thank the reviewer for this suggestion and have implemented it in the text in the style of point 1 of the main reply.

Reply 7:

We have added a sentence below the equation explaining the use of Sklar's theorem in the derivation by the definition of a copula.

Reply 8:

In 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)$.

Reply 9:

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 $\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.

Reply 10:

The 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".

Reply 11:

In the interest of space, we include this as a comment below.

Reply 13:

You 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.

Reply 14:

Please see the pseudo code in the pdf, and our point 2 in the main reply.

Reply 15:

We will include a discussion on the computational cost of our approach. -see reply to other reviewer and experiment on computing time-\ The 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.

Further replies and the proof are in the comment.

Thanks a lot for the detailed clarifications. Trusting the authors to include the proposed changes, I will increase my score to a 5. I am reluctant to increase more, because if I draw a parallel with a journal submission, for the latter I would have liked to review and carefully proofread a revised version of the paper.

We thank the reviewer for their effort in reviewing and offer answers to their comments below.

Reply to weaknesses:

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) 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.

Dimension GMM 1 GMM 2 GMM 3 GMM 4 GMM 5

400

QB-Vine 29.2213 29.7847 29.2775 29.7731 29.0477

RQNSF 29.7247 29.2993 29.7835 29.3515 29.8447

500

QB-Vine 32.7893 33.1789 32.8354 33.4401 32.5520

RQNSF 33.2011 32.6044 33.4355 32.7249 33.4143

600

QB-Vine 35.7948 36.5586 35.8328 36.6095 35.9390

RQNSF 36.4700 35.6756 36.4400 35.7731 36.7090

\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.}

Reply Q1:

Yes, 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.

Reply Q2

We 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.

Reply Q3

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.

Reply Q3

We thank the reviewer for spotting this typo and have fixed it.

Reply Q4

We 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.

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.

Reply Q1:

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.

Reply Q2:

We 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.

Reply Q3:

The 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].

Reply Q4:

We thank you for this comment. We will expand upon this point in the camera ready.

Reply to weakness 1

We 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.

Reply to weakness 2

The 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.

Our 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.

We 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.

Reply to weakness 3:

We 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.

The 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.

Reply Q1:

Thank 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.

Reply Q2:

The 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.