Breaking the Curse with BAND: Nonparametric Distribution Estimation in High Dimensions Shuo-Chieh Huang 1 Chien-Ming Chi 2 Jau-er Chen 3
arXiv:2607.26955v1 [stat.ML] 29 Jul 2026
Abstract
(Liu et al., 2023), distributional random forests (Cevid et al., 2022), and score-based diffusion models (Chen et al., 2023; Benton et al., 2023; Oko et al., 2023; Zhang et al., 2024) achieve minimax-optimal rates but still scale poorly with dimension.
Minimax-optimal rates for multivariate distribution estimation are known to suffer from the curse of dimensionality. We propose a sparse Bayesian network approach in which each conditional probability is estimated using sparsityaware conditional mean methods. The resulting estimator, BAyesian Network Distribution regression (BAND), handles mixed data types in highdimensional time series and achieves polynomial total variation convergence rates while allowing the feature dimension to grow polynomially with the sample size. These rates are substantially faster than the classical optimal rates for multivariate histogram density estimators that lack sparsity. Empirical evaluations show that BAND performs competitively for data sampling and confidence region forecasting against a range of state-of-theart benchmarks.
In contrast, by exploiting the underlying sparsity structure, high-dimensional conditional mean modeling can circumvent such dimensionality issue, using tools such as regularized regression (Tibshirani, 1996), tree-based models (Breiman et al., 2017), and neural networks with bounded fan-in (Lee et al., 1996). Inspired by the success in the aforementioned methods, we aim to tackle the curse of dimensionality in distribution modeling by combining a (sparse) Bayesian network representation and the state-of-the-art conditional mean modeling methods. Using a Bayesian network (Jordan et al., 1999; Bengio & Bengio, 1999), the joint distribution of p variables can be represented as a product of p autoregressive conditional probability functions (see (2) in Section 2.2). This representation underlies many modern architectures in language modeling (Bengio et al., 2003; Child et al., 2019), image generation (Van Den Oord et al., 2016; Salimans et al., 2017), and genetic data analysis (Friedman et al., 2000).
1. Introduction Distribution modeling underlies many statistical tasks, including multivariate confidence region construction (Šidák, 1967), out-of-distribution detection (Liu & Paparrizos, 2024), missing value imputation (Muzellec et al., 2020; Huang et al., 2025), data visualization, and generative modeling. Despite its importance, nonparametric distribution estimation remains difficult in high dimensions. Classical theory shows that, without additional structure, the optimal rate of convergence in total variation is n−β/(2β+p) , where β > 0 is related to the smoothness of the underlying density and p the ambient dimension (Tsybakov, 2009), which is only informative when p = o(log n). This limitation reflects the curse of dimensionality. Recent methods based on kernel estimation (Biroli & Mézard, 2024), adaptive partitioning
In this paper, we propose a sparse Bayesian network model in which each conditional probability is estimated using conditional mean methods that exploit sparse dependence, such as ℓ1 -regularized regression and tree-based models. The resulting estimator, BAyesian Network Distribution regression (BAND), mitigates the curse of dimensionality in highdimensional distribution modeling. The sparse Bayesian network assumption not only allows a formal bias-variance analysis of BAND’s high-dimensional convergence rates, but is also flexible enough to capture complex structures such as mixture distributions. Using one-hot encodings and discretization-based regression methods (Foresi & Peracchi, 1995; Hall et al., 1999; Chernozhukov et al., 2013; Kneib et al., 2023), BAND seamlessly handles discrete, continuous, and mixed data types. In addition, our theoretical results incorporate not only i.i.d. data but also time series data in which case the stationary distribution is of interest. Under a sparse Bayesian network structure, where each variable depends
1
Department of Statistics, Rutgers University, Piscataway, NJ, U.S.A. 2 Institute of Statistical Science, Academia Sinica, Taipei, Taiwan 3 School of Political Science and Economics, National Taiwan University, Taipei, Taiwan. Correspondence to: ChienMing Chi <[email protected]>. Preprint. July 30, 2026.
1
Distribution Modeling via Sparse Bayesian Networks
on at most s0 other variables, BAND achieves polynomial total variation convergence rates with the feature dimension allowed to grow polynomially with the sample size. For continuous time series with strong mixing, the convergence rate is of order p n−1/(3+s0 ) , up to mild subpolynomial factors. This rate is comparable to the state-of-the-art result for high-dimensional density estimation (Vandermeulen et al., 2024). For discrete distributions, we obtain a convergence rate which scales linearly with p under the sparse Bayesian network assumption, which is far faster than the pessimistic minimax convergence rate without sparsity (Han et al., 2015).
1.2. Notation Let (Ω, F, P) be a probability space. Boldface symbols such #– as X denote random vectors, while #– x denotes deterministic vectors. The indicator function is 1{·}, and R is the Borel Pk σ-algebra on R. For #– x ∈ Rk , let ∥ #– x ∥2 = ( j=1 x2j )1/2 and ∥ #– x ∥∞ = maxj |xj |. We denote [1 : p] = {1, . . . , p} #– #– #– #– and write h 1:p = ( h 1 , . . . , h p ). For S ⊂ [1 : p], X S #– denotes the sub-vector of X indexed by S.
2. Mixed-Type Stationary Distribution Modeling #– Let {X t = (X t1 , . . . , X tp )⊤ }nt=1 be the observed pdimensional stationary process (which includes the case of i.i.d. observations). Our goal is to learn its stationary #– distribution. In general, X t may include lagged variables, useful in time series applications (see Section 5). We use #– X = (X 1 , . . . , X p )⊤ to denote a generic random vector #– with joint distribution same as Xt .
We also evaluate BAND empirically on synthetic data, illustrating its performance in data sampling and for forecasting confidence regions in monthly U.S. economic time series. In these experiments, BAND performs competitively with several generative models and density estimators, including normalizing flows (Durkan et al., 2019) and high-dimensional copulas (Nagler & Czado, 2016). Implementations and experiments are available at https: //github.com/no-name213/band.
2.1. Mixed-Type Stationary Distributions 1.1. Related Work Multivariate distribution modeling has accumulated a sizable literature for the fixed-dimensional regime. Copulabased methods often require restrictive assumptions, such as the simplified vine copula (Nagler & Czado, 2016; see also Section 4), or are only supported by simulation studies (Oh & Patton, 2017). Kernel density estimation (Biroli & Mézard, 2024), adaptive partitioning (Liu et al., 2023), and distributional random forests (Cevid et al., 2022) are also classical methods that do not scale to high-dimensional time series data.
Figure 1. (Left) The two-dimensional feature space is partitioned into four regions due to discrete points A and B, including three “discrete” regions indexed by Θ = {({1}, A), ({2}, B), ({1, 2}, (A, B))} and the rest corresponding to the continuous part. (Right) Hypercube partition of the feature support and bins along each coordinate. The discrete points A and B are themselves bins. There are three bins along the first dimension, and, with C as an additional split, also three bins along the second dimension. Together, their Cartesian products constitutes nine hypercubes. Regions beyond the boundary points D and E at (±δn , ±δn ) are omitted.
Recent literature also shed light on when the curse of dimensionality can be avoided using score-based diffusion methods (Song et al., 2020), such as Barron space-valued log-relative-density (Cole & Lu, 2024), low-dimensional support (Oko et al., 2023; Beyler & Bach, 2025), and weak log-concavity (Silveri & Ocello, 2025). Alternatively, Liang et al. (2026) characterizes the complexity of diffusion denoising through a notion of average curvature.
We consider stationary distributions which are possibly of mixed-type. That is, each variable can be continuous, discrete, or mixed (with continuous and discrete components). #– In general, the distribution of X can be defined by several (sub-probability) density functions (given in (1) below). Let Dj denote the set of discrete points for the j-th coordinate, and define Θ = {(S, (dj )j∈S ) S ⊆ [1 : p], dj ∈ Dj for all j ∈ S}. Then, the feature space can be partitioned into subsets indexed by (S, (dj )j∈S ) ∈ Θ, and a subset in which all coordinates are continuous. See the left panel of Figure 1 for an example. To precisely define the densities, we introduce the following decomposition for #– every measurable A ⊆ Rp and each (S, d S ) ∈ Θ,
In comparison, BAND escapes the curse of dimensionality by leveraging the sparse Bayesian network that allows for flexible dependence structures. The sparse Bayesian network is similar to the graph resilience assumption of the recent work of Vandermeulen et al. (2024), which was shown to determine the sample complexity for distribution estimation with i.i.d. data supported on [0, 1]p . However, their work only proves the sample complexity whereas this paper develops an implementable method and is applicable to general time series with mixed-type data.
#– #– #– A(S, d S ) = A ∩ W1 (S, d S ) ∩ W2 (S, d S ), 2
Distribution Modeling via Sparse Bayesian Networks
#– with W1 (S, d S ) = ∩j∈S { #– v ∈ Rp : vj = dj }, #– #– W2 (S, d S ) = ∩j∈S c { v ∈ Rp : vj ̸∈ Dj }, S c := [1 : #– p] \ S, and Ac = A \ ∪(S, #– d S )∈Θ A(S, d S ) . Then the #– #– stationary distribution of X is represented by P(X ∈ A) = P #– #– #– P(X ∈ Ac ) + (S, #– d S )∈Θ P(X ∈ A(S, d S )), where Z #– P(X ∈ Ac ) = πc ( #– x ) d #– x,
#– where ωj : {0, 1}mj → {1, . . . , mj } satisfies ωj ( h j ) = q #– if the qth coordinate of h j is 1. For notational complete#– # – #– #– #– #– ness, we write P(H 1 = h 1 | H 0 = h 0 ) = P(H 1 = h 1 ) #– and e1l ( h 0 ) = E(H 1l ). Now, we can estimate each ejl using high-dimensional regression techniques such as decision trees (Morgan & Sonquist, 1963; see Section 3.2), CART (Breiman et al., 2017), or Lasso (Tibshirani, 1996). Let ebjl denote the estimated probability function correspond#– ing to ejl . Then the estimated probability P(X ∈ A) for a measurable set A ⊂ Rp is then defined as
Ac
(1) Z #– #– P(X ∈ A(S, d S )) #– #– #– ϕS, d S ( u ) d u , = #– #– #– #– P(X S = d S ) u ∈Ā(S, d S ) #– #– in which Ā(S, d S ) = { #– a S c : #– a ∈ A(S, d S )}, and let c #S πc : Rp → [0, ∞) and ϕS, #– → [0, ∞) denote dS : R the corresponding conditional density functions. For notational completeness, the second integral in (1) equals 1 when S c = ∅, in which case ϕS, #– d S ≡ 1. We further define the p extended conditional density πS, #– d S : R → [0, ∞) such #– #– #– #– #– that πS, d S ( x ) = ϕS, d S ( x S c ) for all x ∈ Rp , rendering it invariant to the discrete coordinates indexed by S.
µ b(A) =
#– h 1:p
ν #– h 1:p (A)
p Y
#– ebj,ωj (hj ) ( h j−1:1 ),
(3)
j=1
#– where the sum is over all one-hot vectors h j ∈ {0, 1}mj , and that ν #– h 1:p is a normalized Lebesgue measure with re#– spect to the continuous coordinates of Cube( h 1:p ), the hypercube corresponding to the ωj (hj )-th bin on the j-th coordinate, which has total mass one and zero outside the hypercube. For example, for the hypercube whose second coordinate is fixed at B in Figure 1, its continuous coordi#– nate is the first coordinate. In addition, if Cube( h 1:p ) is #– a discrete cube, then ν #– h 1:p (A) = 1 if Cube( h 1:p ) ⊆ A and ν #– (A) = 0 otherwise. Since this method relies on h 1:p the Bayesian network representation as well as modeling multiple probability functions via regressions, we call it the BAyesian Network Distribution (BAND) regression.
2.2. Bayesian Network Distribution Regression We employ a discretization strategy to learn the mixedtype, multivariate stationary distribution. First, each coordinate is partitioned into mj bins and the one-hot vectors #– H j = (H j1 , . . . , H jmj )⊤ encodes which bin Xj lies in. Note that the mj th bin is reserved to represent the tails Ij \ [−δn , δn ], where Ij denotes the support of X j with P(X j ∈ Ij ) = 1, and the remaining (mj − 1) bins form a partition of Ij ∩ [−δn , δn ]. In addition, we also treat the discrete points in Dj as independent bins. Then, the Cartesian products of the bins in each coordinate constitute the resulting p-dimensional hypercubes, which form a partition of the feature support. An example is illustrated in the right panel of Figure 1. These hypercubes are used in our method (to be introduced below) through the one-hot encodings. For now, we treat them as given. The optimal bin sizes are discussed in Section 3, and we also tune them using some practical methods in Sections 4 and 5.
Many existing methods for discrete distribution modeling use Bayesian networks via maximizing the log-likelihood function, using shallow neural networks such as sigmoid belief networks (Neal, 1992; Jordan et al., 1999), or deep architectures such as PixelRNN (Van Den Oord et al., 2016) and sparse transformers (Child et al., 2019). In contrast, the proposed BAND regression models the distribution by directly estimating each conditional probability function #– #– #– #– ejl (H j−1 , . . . , H 1 ) = E H jl | H j−1 , . . . , H 1 . A key advantage of this formulation is that BAND explicitly handles variable sparsity structure using the well-known highdimensional regression techniques, which provides provable guarantees that avoid the curse of dimensionality.
We now introduce the Bayesian network representation (Jordan et al., 1999; Bengio & Bengio, 1999) which is fundamental to our method. Specifically, for every Pp #– #– #– ( h 1 , . . . , h p ) ∈ {0, 1} l=1 ml , where each h l is a one-hot vector, it holds that P
X
3. Theoretical Foundations of Distribution Modeling
p \ #– #– {H j = h j }
3.1. Convergence Rates in Total Variation Distance
j=1
=
=:
p Y j=1 p Y
#– #– P Hj = h j
j−1 \
#– #– {H l = h l }
Our main goal in this subsection is to derive a convergence rate in total variation distance (Theorem 3.3) for BAND. To this end, let us first introduce some high-level assumptions. #– Condition 3.1. There is dmin > 0 such that ejl ( h j−1:1 ) ≥ #– #– dmin if ω1 ( h 1 ) < m1 , . . . , ωj−1 ( h j−1 ) < mj−1 and l < mj . Additionally, P(E † ) ≤ ∆ for
(2)
l=1
#– #– ej,ωj ( #– h j ) ( h j−1 , . . . , h 1 ),
j=1
3
Distribution Modeling via Sparse Bayesian Networks
some ∆ > 0, where E † = ∩1≤j≤p ∩1≤l<mj √ #– #– ∩ωq ( #– bjl ( h j−1:1 ) ≤ T . h q )<mq ,1≤q<j ejl ( h j−1:1 ) − e
polynomial convergence rates in n, allowing the feature dimension p to grow polynomially with the sample size. In contrast, existing theoretical results achieving optimal rates (Cevid et al., 2022; Liu et al., 2023; Oko et al., 2023; Zhang et al., 2024; Biroli & Mézard, 2024) typically require p ≤ log n for total variation consistency.
√ Condition 3.1 assumes a deviation upper bound T for the regression estimates, and a probability lower bound dmin for the regression functions. These parameters can be estimated under more fundamental assumptions (see Corollaries 3.7 √ and 3.9), with T depending only mildly on p thanks to the sparse Bayesian networks (Condition 3.5). #– In Condition 3.2 below, Θ⋆ = Θ if X is discrete; otherwise, Θ⋆ = Θ ∪ {c} includes an additional label c. Recall that πc is defined in (1) while πθ , θ ∈ Θ, are defined after (1).
3.2. Application to Continuous Multivariate Time series #– Suppose {X t }nt=1 stationary process with R is a continuous Θ = ∅ (that is, Rp πc ( #– x ) d #– x = 1) and feature supports #– I1 = · · · = Ip = R. Let H tj = (H tj1 , . . . , H tjmj )⊤ , #– t ∈ [1 : n], denote the sample counterparts of H j , i.e., the corresponding discretized sample constructed as in Section 2.2. To illustrate how BAND adapts to the sparse Bayesian network, consider, for each 1 ≤ j ≤ p and 1 ≤ l < mj , a sparse tree model
Condition 3.2. There exists some R > 0 such that for each θ ∈ Θ⋆ and every #– x ∥∞ ≤ R, it holds that x with ∥ #– πθ ( #– x ) > 0, log πθ ( #– x ) ∈ C 1 ( #– x ), and ∥∇ log πθ ( #– x )∥2 ≤ √ K1 p R for some constant K1 > 0.
#– #– ebjl ( h j−1 , . . . , h 1 ) =
A bounded gradient of the log-density (the score) is a standard regularity assumption in the density estimation literature, including in score-based diffusion models. Condition 3.2 additionally accounts for the dependence on the feature dimension p. It can be verified that the (boundary adjusted) uniform, Student-t, and multivariate Gaussian distributions and their mixtures all satisfy Condition 3.2; see Section D.1 for details. Alternatively, one may assume ∥∇ log πθ ( #– x )∥∞ ≤ K1 R, which implies Condition 3.2 √ since ∥∇ log πθ ( #– x )∥2 ≤ p∥∇ log πθ ( #– x )∥∞ .
Pn
#–
#–
t=1 H tjl
Q
Pn
q∈Sj 1{H tq = h q }
t=1
Q
q∈Sj 1{H tq = h q }
#–
#–
(4)
#– #– for ( h 1 , . . . , h p ) ∈ {0, 1}m1 ×· · ·×{0, 1}mp with one-hot #– #– h j and ωj ( h j ) < mj ; otherwise, the tree model prediction is defined to be zero. Here, S1 , . . . , Sp are relevant sets given by Condition 3.5 below. Condition 3.5 (Sparse Bayesian Networks). (a) For each #– j, there is some Sj ⊂ {1, . . . , j − 1} such that H j is inde#– #– pendent of (H l , l ̸∈ Sj , l < j) conditional on (H l , l ∈ Sj ). (b) There are S1 , . . . , Sp such that Sj ⊂ {1, . . . , j − 1} and #– #– that X Sj ∪{j} is independent of X {1,...,j−1}\Sj .
Let {δn } be a positive sequence, and let εn ≥ 0 denote the maximum diameter of the cubes in the partition of [−δn , δn ]p (cubes as defined after (3)). Now we are ready to state the main result, Theorem 3.3, whose proof is provided in Section B.1. Note that (∆, dmin , T, E † , ε, p) all depend on n; however, this dependence is omitted for brevity. #– #– #– Theorem 3.3. Set ebjl ( h j−1 , . . . , h 1 ) = 0 if ωq ( h q ) = mq for some 1 ≤ q < j or l = mj . For all large n, assume Condition 3.1–3.2 with R = δn . Assume also ε ≤ √ √ (δn pK1 )−1 and ep T ≤ dmin for all large n. Then, it holds for all large n that, with probability at least 1 − ∆,
Condition 3.5 is trivially satisfied when the Bayesian network is not sparse (Sj = {j − 1, . . . , 1}). Condition 3.5(a) requires the discretized one-hot vectors to satisfy an autoregressive conditional independence property, which does not automatically follow from X 1 , . . . , X p having the same property. However, it is a natural condition in many applied fields, such as economics, finance, and engineering, where continuous variables are routinely discretized. For example, Tauchen (1986) shows that increasingly fine partitions yield accurate approximations of continuous processes by finite-state Markov chains. Additionally, Example 4.2 in Section 4 illustrates that, with sufficiently many discretization splits, Condition 3.5(a) holds for some mixtures of uniform distributions. More examples can be found in (Vandermeulen et al., 2024). Alternatively, Condition 3.5(b) assumes there are multiple groups of mutually independent variables, which is another form of sparse Bayesian networks.
√ √ #– ep T sup P(X ∈ A) − µ b(A) ≤ + ε pδn K1 e + Qδn , p d min A∈R
#– where Qδn = P(X ̸∈ [−δn , δn ]p ). Remark 3.4. The first two terms in the upper bound in Theorem 3.3 correspond to estimation variance and bias, respectively. In Section B.1, we also provide a similar result for conditional distribution. One major implication of Theorem 3.3 is that BAND can consistently estimate high-dimensional distributions, which allows p to grow polynomially with the sample size n. To illustrate, in the following subsections we consider the special cases of continuous and discrete distributions in Corollaries 3.7 and 3.9, respectively. Across both cases, we achieve
With Condition 3.5, BAND combined with the sparse tree models (4) can potentially attain a faster rate of convergence; see Corollaries 3.7 and 3.9, which are derived assuming these relevant sets are known. A formal analysis of recovering these sets via sparse learning is left for future 4
Distribution Modeling via Sparse Bayesian Networks
work. In practice, the relevant sets Sj are unknown and can be estimated using sparse regression methods such as Lasso (Tibshirani, 1996), CART (Breiman et al., 2017), and CatBoost (Prokhorenkova et al., 2018).
ate histogram estimators (Beirlant & Gyorfi, 1998). By comparison, the rate derived in Example 3.8 circumvents the curse of dimensionality by replacing the dimension p in the exponent with the sparsity-related term 1 + s0 ≤ p (where s0 = 0 corresponds to the case Sj = ∅ for all j), with an additional linear factor of p. Vandermeulen et al. (2024) obtained convergence rates comparable to that in Example 3.8, specifically by replacing the dimension p in the exponent with a sparsity parameter under a condition similar to Condition 3.5. However, their results depend on densities supported exclusively on the unit hypercube [0, 1]p . Such confined support is generally incompatible with time series data, which typically requires a broader domain, thus limiting their findings to independent data settings.
Finally, we need a technical, though standard, assumption on the weak dependence in the data in order to establish concentration inequalities. Condition 3.6 (Strong Mixing). There are constants γ0 , γ1 > 0 such that the sequence {(X tl , l ∈ Sj ∪ {j})}t is strongly mixing with mixing coefficients satisfying α(k) ≤ γ0 e−γ1 k for all k ≥ 1 and all 1 ≤ j ≤ p, with Sj specified in Condition 3.5. A formal definition of the strong mixing property is provided in Section D.3. Condition 3.6 requires that each sub-process {(X tl , l ∈ Sj ∪{j})}t is strong mixing, allowing linear AR, threshold AR, bounded nonlinear AR, ARCH, GARCH, and VAR processes (An & Huang, 1996; Tsay, 2005). This is slightly weaker than requiring the full process to be strongly mixing. Now we can state Corollary 3.7, whose proof is given in Section B.2.
Besides sparse Bayesian networks, prior literature on copula models (Nagler & Czado, 2016; Oh & Patton, 2017) have considered other dependence structures for multivariate distribution modeling, such as the simplified vine copulas. However, those assumptions are rather restrictive, which, for example, exclude some multi-modal mixture distributions. In contrast, the sparse Bayesian network condition employed here has greater flexibility.
Corollary 3.7. Let δn > 0 be arbitrary. Assume Condition 3.2 (with R = δn ) and Condition 3.6, and either Condition 3.5(a) or 3.5(b) with max1≤j≤p #Sj ≤ s0 for some constant s0 ≥ 0. In addition, suppose the joint #– density of X J is uniformly bounded above and below whenever #J ≤ s0 + 1. If the splitting scheme of the bins satisfies max1≤j≤p mj ≤ n, supn≥1 (L2n /L1n ) < 0 ∞, L−1−s ≤ n, and limn→∞ L2n = 0, where 1n L1n and L2n denote the minimum and maximum bin s0 1 −1/2 5 − 2 − 2 −b0 lengths, and if pn (log n) L2n ≤ 1 and pL2n δn = o(1) for some b0 > 0, then with probabil#– ity tending to one, supA∈Rp P(X ∈ A) − µ b(A) ≤ 1
s
− 1 − 20 −b0
pn− 2 (log n)5 L2n2
3.3. Application to High-Dimensional Discrete Data #– In this subsection, we consider a discrete process X t whose stationary distribution satisfies P(X j ∈ Dj ) = 1 with Dj = {1, . . . , mj −1}, and define Ij := Dj ∪{m0 }, where m0 = max1≤j≤p mj and δn = m0 − 1 are assumed to be constant. Then, by the construction of the bins introduced in Section 2.2, X j = l if and only if H jl = 1 for l ∈ [1 : (mj − 1)]. Here, all cubes here are discrete cubes, with diameter ε = 0. Corollary 3.9 builds on the same tree model (4); see Section B.3 for the proof.
+ pL2n δn log n + Qδn . Corollary 3.9. Let constants s0 ≥ 0 and γ3 , γ4 , γ5 > 0 be given. Assume Condition 3.5(a) with max1≤j≤p #Sj ≤ s0 , Condition 3.6, and that max1≤j≤p mj ≤ m0 and −1/2 pnQ (log n)5 ≤ 1. In addition, assume γ4 γ3#S ≤ E( j∈S H jlj ) ≤ γ5 γ3#S for all 1 ≤ lj < mj and S ⊂ {1, . . . , p} with #S ≤ s0 +1. Then, with probability tending #– to one, supA∈Rp P(X ∈ A) − µ b(A) ≤ pn−1/2 (log n)5 .
Note that Corollary 3.7 allows the feature dimension p to grow polynomially with the sample size, while maintaining polynomial convergence rates. Leveraging the sparsity in Condition 3.5, we avoid the curse of dimensionality and achieve a convergence rate far faster than the existing results for adaptive partitioning (Liu et al., 2023), distributional random forests (Cevid et al., 2022), and score-based diffusion models (Oko et al., 2023; Zhang et al., 2024), when p ≫ log n. In Example 3.8 below, we explicitly calculate the polynomial convergence rate with a carefully selected bin size.
In general, the minimax optimal rate pfor estimating discrete distributions without sparsity is O( Support Size/n) (Han et al., 2015). Take a p-dimensional vector of independent Bernoulli variables as an example. The support size is 2p , which scales poorly with p (in fact, the minimax rate is much slower when the support size increases with n). In comparison, Corollary 3.9 offers a substantially faster convergence rate that scales linearly with p by incorporating the sparse Bayesian network structure.
Example 3.8. Let L2n = n−1/(3+s0 +2b0 ) and δn = nϵ1 for some ϵ1 > 0, assuming πc satisfies sufficiently strong polynomial moment conditions so that Qδn is negligible. Then the convergence rate in Corollary 3.7 is p n−1/(3+s0 +2b0 )+ϵ1 (log n)5 . The rate n−1/(2+p) is the established optimal for multivari5
Distribution Modeling via Sparse Bayesian Networks
4. Simulation Experiments
the jth variable from a hypercube determined by an mj dimensional probability vector #– p conditional on the preceding j − 1 variables (see Section A.2). To reduce instability, we introduce a sampling probability threshold factor (SPTF ∈ [0, 0.4]), zeroing entries of #– p below SPTF × std( #– p ); similar techniques have been used in generative models for sampling stabilization (Lu et al., 2025).
4.1. Simulation Setup In this section, we first examine the finite-sample performance of BAND when applied to i.i.d. synthetic data generated from bimodal distributions described in Examples 4.1– 4.2 below. Multimodal distributions arise in many applications, from discrete mixtures in genetics (Friedman et al., 2000) and text classification (Juan & Vidal, 2002) to continuous regime-switching models (Tauchen, 1986; Hamilton, 1989) and clustering (Melnykov & Maitra, 2010), but learning multimodal laws is difficult in high dimensions. In the following, let #– µ 1 = (0, . . . , 0)⊤ and #– µ 2 = (r, . . . , r)⊤ , where r ≥ 0. #– Example 4.1 (Gaussian mixture). X ∼ P2 1 #– i=1 2 Np ( µ i , Ip ). #– ∼ Example 4.2 (Uniform mixture). Let X P2 1 #– #– #– #– 3 3 p U ([− , ] + µ ), where A + v = { x + v | p i i=1 2 2 2 #– x ∈ A} for A ⊂ Rp and #– v ∈ Rp .
Since we use CatBoost as the sub-routine for BAND, the computational cost of training BAND is proportional to training (p n1/3 )2 CatBoost models on an n × p feature matrix; see Prokhorenkova et al. (2018) for CatBoost’s complexity. In typical economic time-series applications, the number of variables is modest (see Section 5), so the computational cost is manageable. Moreover, we found that the default learning rates in the range [0.03, 0.05] for CatBoost perform as stably as automatic tuning. Therefore, we simply fix the learning rate at 0.05 in the experiments. Finally, note that tuning SPTF does not require retraining the model, which is computationally inexpensive. Gaussian mixture models (Scrucca et al., 2016): Gaussian mixtures models are fitted using the mclust package in R, which implements the expectation–maximization algorithm with the number of components fixed at K = 2.
Note that Examples 4.1–4.2 satisfy the sparsity Condition 3.5(a) with Sj = ∅ for r = 0. Meanwhile, for Example 4.2 with r > 3, Condition 3.5(a) is met with Sj = {j − 1}, provided there are sufficiently many evenlyspaced splits and δn ≥ 3/2 + r, resulting in two wellseparated modes (see Figure 2); see Section D.2 for details. In general, however, the sparsity condition does not hold, creating challenges for BAND in learning the distribution.
Vine Copula (Nagler & Czado, 2016): We implement the vine copula model via the kdvine package in R, which was proposed to estimate multivariate densities. Normalizing Flow (Durkan et al., 2019): Normalizing flows are mature and well-validated generative models, which implicitly learns the multivariate density. Compared to the data-hungry score-based diffusion models (Song et al., 2021; Zhang et al., 2025), it offers a practical and strong benchmark. We implement normalizing flows using the nflows library (Durkan et al., 2020). To maintain robustness in our moderate-sample setting, we vary only the number of flow layers and use a standard 70/30 training-validation split for model selection, keeping all other hyperparameters at the library’s stable defaults.
Next, we consider the following block-wise graphical Gaussian distribution. Example 4.3. Partition {1, . . . , p} into m equal-sized #– #– blocks B1 , . . . , Bm . The data follow X Bk ∼ N mp ( 0 , (1 − #– #– ρ)I mp + ρ 1 1 ⊤ ), drawn independently for k = 1, 2, . . . , m, #– where 1 denotes a vector of ones. Note that Condition 3.5(b) holds in Example 4.3, where a larger m indicates a higher degree of sparsity. The sparsity pattern is not known to the practitioner. Nevertheless, as we will see later, BAND can successfully adapt to the sparse structure and estimate the underlying distribution well.
4.2. Results To assess the performance of the employed methods, we use the Cramér–Wold (CW) distance (Cramér & Wold, 1936), a computable lower-bound surrogate for the total variation distance; details are provided in Section A.1. Tables 1–3 show the resulting CW distances for various (p, n, r) configurations.
4.1.1. I MPLEMENTATION AND B ENCHMARK M ETHODS Here we briefly discuss the implementation details and introduce some benchmark methods. BAND: For continuous variables, as in our applications, we recommend using ⌈factor × n1/3 ⌉ equally-spaced breaks over the sample range of each coordinate, following the classical cube-root rule for binning (Scott, 1979). We use factor = 2 in this section. The mean regression models are fitted using CatBoost (Prokhorenkova et al., 2018) with default parameters due to its minimal hyperparameter tuning (McElfresh et al., 2023). When sampling, BAND draws
For Example 4.1, mclust is correctly specified and therefore serves as the ideal benchmark method. As shown in Table 1, mclust yields CW errors that are among the lowest. Nevertheless, BAND fares quite well, which significantly outperforms nflows in the high dimension case p = 30.
6
Distribution Modeling via Sparse Bayesian Networks
For Example 4.2 in Table 2, BAND yields the lowest CW errors when n = 1000. In particular, when the two supports in the uniform mixture are well separated (r = 4), the CW error of BAND (around 0.055) is 75% smaller compared to kdvine (around 0.20–0.25). Indeed, Figure 2 shows that BAND accurately recovers the uniform mixture distribution, while kdvine is unsuitable for this case. nflows also incurs high CW errors despite its flexibility in distribution estimation. We note that although mclust is misspecified in this example, its CW errors in Table 2 do not appear particularly poor. Table 1. Average Cramér–Wold distance over 30 trials for the four methods. Data are generated from the Gaussian mixture in Example 4.1. Standard deviations are in parentheses. Method
p/n
r=0
r=2
r=4
BAND
10/100 10/1000 30/100 30/1000
0.169 (0.02) 0.055 (0.01) 0.191 (0.02) 0.062 (0.01)
0.158 (0.02) 0.070 (0.01) 0.180 (0.02) 0.077 (0.01)
0.161 (0.03) 0.059 (0.01) 0.181 (0.02) 0.062 (0.01)
mclust
10/100 10/1000 30/100 30/1000
0.131 (0.03) 0.057 (0.01) 0.153 (0.02) 0.067 (0.01)
0.117 (0.02) 0.055 (0.01) 0.151 (0.03) 0.065 (0.01)
0.134 (0.02) 0.057 (0.01) 0.147 (0.02) 0.063 (0.01)
nflow
10/100 10/1000 30/100 30/1000
0.202 (0.04) 0.087 (0.01) 0.409 (0.05) 0.087 (0.01)
0.182 (0.02) 0.090 (0.01) 0.358 (0.05) 0.093 (0.01)
0.187 (0.02) 0.106 (0.02) 0.280 (0.04) 0.122 (0.02)
kdvine
10/100 10/1000 30/100 30/1000
0.166 (0.02) 0.076 (0.01) 0.176 (0.02) 0.090 (0.01)
0.174 (0.03) 0.150 (0.01) 0.202 (0.02) 0.183 (0.01)
0.219 (0.03) 0.198 (0.01) 0.260 (0.02) 0.244 (0.01)
p/n
r=0
r=2
r=4
BAND
10/100 10/1000 30/100 30/1000
0.182 (0.02) 0.059 (0.01) 0.195 (0.02) 0.065 (0.01)
0.157 (0.02) 0.063 (0.01) 0.184 (0.02) 0.069 (0.01)
0.147 (0.02) 0.054 (0.01) 0.164 (0.02) 0.059 (0.01)
mclust
10/100 10/1000 30/100 30/1000
0.137 (0.02) 0.088 (0.01) 0.148 (0.01) 0.092 (0.00)
0.114 (0.02) 0.072 (0.01) 0.147 (0.02) 0.076 (0.01)
0.114 (0.02) 0.061 (0.00) 0.139 (0.02) 0.065 (0.01)
nflow
10/100 10/1000 30/100 30/1000
0.188 (0.02) 0.103 (0.02) 0.375 (0.05) 0.116 (0.01)
0.181 (0.03) 0.088 (0.01) 0.318 (0.04) 0.102 (0.01)
0.185 (0.03) 0.113 (0.02) 0.254 (0.05) 0.124 (0.02)
kdvine
10/100 10/1000 30/100 30/1000
0.167 (0.03) 0.079 (0.01) 0.175 (0.02) 0.089 (0.01)
0.167 (0.01) 0.161 (0.01) 0.195 (0.01) 0.190 (0.01)
0.212 (0.01) 0.206 (0.01) 0.257 (0.01) 0.251 (0.01)
(b) mcluster
(c) nflow
(d) kdvine
Figure 2. Scatterplots of 10,000 empirical data points (blue) and 10,000 points sampled from the four fitted models (red). Data are drawn from Example 4.2 with r = 4. Each method is trained on n = 1000 samples with p = 2. Table 3. Average Cramér–Wold distance over 30 trials for the four methods. Data are generated from as in Example 4.3. We set p = 30, ρ = 1/2, n = 1,000. Standard deviations are in parentheses. #blocks (m) BAND mcluster nflow kdvine
3
6
10
30
0.082 (0.009) 0.065 (0.007) 0.091 (0.011) 0.084 (0.007)
0.084 (0.007) 0.064 (0.007) 0.094 (0.010) 0.088 (0.009)
0.077 (0.008) 0.065 (0.008) 0.095 (0.013) 0.087 (0.009)
0.061 (0.007) 0.065 (0.007) 0.095 (0.015) 0.088 (0.007)
5. Application: Multivariate Forecasting Confidence Regions for Macroeconomic Time Series
Table 2. Average Cramér–Wold distance over 30 trials for the four methods. Data are generated from the uniform mixture in Example 4.2. Standard deviations are in parentheses. Method
(a) BAND
In this section, we apply BAND to construct (multivariate) forecasting confidence regions for U.S. macroeconomic time series. Macroeconomic forecasting is important but challenging, as economic time series often exhibit nonstationarity, nonlinear dynamics, and may contain substantial outliers. We focus on three key monthly series, namely unemployment rate, federal funds rate, and inflation rate. Before analysis, it is standard practice to remove unit root nonstationarity (McCracken & Ng, 2016). We analyze URt , the first difference of the monthly U.S. unemployment rate, FRt , the first difference of the monthly effective federal funds rate, and IRt , the monthly inflation rate computed as IRt = (CPIt −CPIt−1 )/CPIt−1 , where CPIt is the U.S. consumer price index. We consider data from February 1959 to August 2025, which are publicly available from the FRED database (McCracken & Ng, 2016).
Finally, Table 3 reports the CW errors for the graphical Gaussian distribution in Example 4.3. We again use the correctly-specified mclust as the ideal benchmark. Note that BAND successfully adapts to the underlying sparsity and achieves lower CW errors when there are many independent blocks among the p = 30 variables. In contrast, the CW errors of kdvine and nflows remain the same regardless of the sparsity level.
We consider two forecasting periods, Jan. 2014–Dec. 2019 (test period I) and Aug. 2022–Aug. 2025 (test period II) using all observations prior to each period as the corresponding training sample. We train BAND regression models on data #– #– #– #– X t = (U t , U t−k ), where U t = (FRt , IRt , URt ) and the forecasting horizon k ∈ {1, 3, 6}. Table 4 shows some summary statistics of the data. Note that the training data for 7
Distribution Modeling via Sparse Bayesian Networks Table 4. Summary statistics of the training data in the two forecasting periods, including the ranges, interquantile ranges (IQR), bin size (BS), and the BS-to-IQR (BS / IQR) ratio for each series. The selected factor and the resulting number of bins (#Bins) for each forecasting period are also reported.
test period II contain significant outliers in the unemployment series due to COVID-19, as shown by the extremely large range and the inflated bin-size-to-interquantile-range ratio. The outliers could cause instability in conditional distribution estimation. Hence we discuss test period I and II separately in our analysis.
Test Period I Test Period II Jan. 2014–Dec. 2019 Aug. 2022–Aug. 2025 URt FRt IRt URt FRt IRt
The k-month-ahead forecasting confidence regions of BAND at α level are constructed as follows. Suppose #– U t−k = #– u t−k is observed. Then, focusing on the hypercubes whose last three coordinates contain #– u t−k , we select those with the highest predicted probabilities according to BAND until the cumulative probability just exceeds α. Next, as is standard for partition-based methods (Hyndman, 1996), the size of the final selected hypercube is adjusted by proportionally shrinking each dimension toward its center so that the resulting cumulative probability matches α exactly. When constructing confidence regions for univariate variables (i.e., we are interested in marginal confidence regions), coordinates that are not of direct interest are “integrated out” by aggregation. Further details are provided in Section A.3 in the supplementary material. Finally, the bin-size parameter factor (see Section 4) is selected over a grid {0.4, 0.5, . . . , 2.0} by minimizing the maximum discrepancy between the empirical coverage rate and the target α = 0.9 level for the 6-month-ahead forecast #– of U t . Selected factor for each test period is reported in Table 4.
Range 9.69 0.0358 IQR 0.24 0.0033 BS 1.38 0.0051 BS / IQR 5.77 1.55 #Bins 7 factor 0.9
1.60 0.20 0.23 1.15
9.69 0.0358 0.20 0.0033 1.07 0.0040 5.35 1.21 9 1.0
12.60 0.20 1.40 7.00
Table 5. Average absolute gaps (Ave Gap), across different forecasting horizons k ∈ {1, 3, 6}, between the empirical coverage rate and the specified confidence level α. (A) average over univariate forecasts; (B) multivariate forecasts. Test Period
Ave Gap α ≤ 0.9
I α = 0.9
α ≤ 0.9 II
As benchmarks, we employ the quantile regression (QR) and the quantile random forests (QRF, Meinshausen & Ridgeway, 2006) in our analysis. To construct confidence re1+α gions for QR and QRF, we use the two-sided ( 1−α 2 , 2 ) quantile intervals, which are the default in quantreg and quantregForest in R. Multi-dimensional confidence regions are constructed by intersecting the corresponding marginal intervals, each adjusted via the Bonferroni correction: αadj = 1− Region1−α Dimension . Although the Bonferroni correction is known to be conservative and may be less suited for constructing multivariate forecasting confidence regions, such issue has received relatively limited attention in both classical and recent deep-learning approaches (Kollovieh et al., 2023; Kotelnikov et al., 2023).
α = 0.9
Method
(A)
(B)
BAND QR QRF BAND QR QRF
0.114 0.118 0.114 0.060 0.059 0.061
0.066 0.181 0.143 0.032 0.061 0.035
BAND QR QRF BAND QR QRF
0.186 0.115 0.100 0.057 0.082 0.047
0.085 0.266 0.180 0.048 0.089 0.044
suggests that BAND is effective for multivariate distribution estimation when trained on the full set of variables, leveraging its ability to handle higher-dimensional data. For test period II, BAND shows larger average gaps in forecasting individual series because of the outliers. For multivariate forecast confidence regions, BAND and QRF perform comparably at α = 0.9, while BAND achieves approximately 50% smaller average gaps than QRF when averaged over α ≤ 0.9. Both methods outperform QR in these cases. BAND performs more stably in multivariate than in univariate settings in period II. Intuitively, including more variables distributes the probability mass across multiple hypercubes, reducing the chance that a single hypercube dominates, which is a common issue in univariate cases and in the presence of outliers.
5.1. Results We examine the gaps between the realized coverage rates and the target values of α for α ∈ {0.1, 0.2, . . . , 0.9}. Table 5 reports the average gaps for univariate and multivariate forecasts. For test period I, when used to construct marginal forecasting confidence regions, the average coverage gaps of BAND for the univariate series are similar to those of QR and QRF, although it is trained using all six variables. For multivariate forecasting confidence regions (column (B) in Table 5), BAND clearly outperforms QR and QRF. This
In summary, BAND demonstrates strong empirical performance for multivariate distribution estimation and is applicable to mixed-type distributions. By adapting to underlying sparsity, it achieves fast convergence rates that mitigate the curse of dimensionality, making it effective in highdimensional distribution modeling. 8
Distribution Modeling via Sparse Bayesian Networks
Acknowledgements
Chernozhukov, V., Fernández-Val, I., and Melly, B. Inference on counterfactual distributions. Econometrica, 81 (6):2205–2268, 2013.
We thank Tengyuan Liang and Ruey Tsay for helpful comments and constructive suggestions.
Child, R., Gray, S., Radford, A., and Sutskever, I. Generating long sequences with sparse transformers. arXiv preprint arXiv:1904.10509, 2019.
Impact Statement This paper presents work whose goal is to advance the field of Machine Learning. There are many potential societal consequences of our work, none which we feel must be specifically highlighted here.
Cole, F. and Lu, Y. Score-based generative models break the curse of dimensionality in learning a family of subgaussian distributions. In International Conference on Learning Representations (ICLR), 2024.
References
Cramér, H. and Wold, H. Some theorems on distribution functions. Journal of the London Mathematical Society, 1(4):290–294, 1936.
An, H. and Huang, F. The geometrical ergodicity of nonlinear autoregressive models. Statistica Sinica, pp. 943–956, 1996.
Durkan, C., Bekasov, A., Murray, I., and Papamakarios, G. Neural spline flows. In Advances in Neural Information Processing Systems (NeurIPS), 2019.
Beirlant, J. and Gyorfi, L. On the lrerror in histogram density estimation: The multidimensional case. Journal of Nonparametric Statistics, 9(2):197–216, 1998.
Durkan, C., Bekasov, A., Murray, I., and Papamakarios, G. nflows: normalizing flows in PyTorch, November 2020. URL https://doi.org/10.5281/ zenodo.4296287.
Bengio, Y. and Bengio, S. Modeling high-dimensional discrete data with multi-layer neural networks. Advances in Neural Information Processing Systems (NeurIPS), 12, 1999.
Foresi, S. and Peracchi, F. The conditional distribution of excess returns: An empirical analysis. Journal of the American Statistical Association, 90(430):451–466, 1995.
Bengio, Y., Ducharme, R., Vincent, P., and Jauvin, C. A neural probabilistic language model. Journal of Machine Learning Research, 3:1137–1155, 2003.
Friedman, N., Linial, M., Nachman, I., and Pe’er, D. Using Bayesian networks to analyze expression data. In The Annual International Conference on Computational Molecular Biology, pp. 127–135, 2000.
Benton, J., De Bortoli, V., Doucet, A., and Deligiannidis, G. Nearly d-linear convergence bounds for diffusion models via stochastic localization. arXiv preprint arXiv:2308.03686, 2023.
Hall, P., Wolff, R. C., and Yao, Q. Methods for estimating a conditional distribution function. Journal of the American Statistical Association, 94(445):154–163, 1999.
Beyler, E. and Bach, F. Optimal denoising in score-based generative models: The role of data regularity. arXiv preprint arXiv:2503.12966, 2025.
Hamilton, J. D. A new approach to the economic analysis of nonstationary time series and the business cycle. Econometrica, pp. 357–384, 1989.
Biroli, G. and Mézard, M. Kernel density estimators in large dimensions. arXiv preprint arXiv:2408.05807, 2024. Breiman, L., Friedman, J., Olshen, R. A., and Stone, C. J. Classification and regression trees. Chapman and Hall/CRC, 2017.
Han, Y., Jiao, J., and Weissman, T. Minimax estimation of discrete distributions under ℓ1 loss. IEEE Transactions on Information Theory, 61(11):6343–6354, 2015. doi: 10.1109/TIT.2015.2478816.
Cevid, D., Michel, L., Näf, J., Bühlmann, P., and Meinshausen, N. Distributional random forests: Heterogeneity adjustment and multivariate distributional regression. Journal of Machine Learning Research, 23(333):1–79, 2022.
Huang, S.-C., Liang, T., and Tsay, R. S. Temporal wasserstein imputation: Versatile missing data imputation for time series. arXiv preprint arXiv:2411.02811, 2025. Hyndman, R. J. Computing and graphing highest density regions. The American Statistician, 50(2):120–126, 1996.
Chen, S., Chewi, S., Li, J., Li, Y., Salim, A., and Zhang, A. R. Sampling is as easy as learning the score: Theory for diffusion models with minimal data assumptions. International Conference on Learning Representations (ICLR), 2023.
Jordan, M. I., Ghahramani, Z., Jaakkola, T. S., and Saul, L. K. An introduction to variational methods for graphical models. Machine Learning, 37(2):183–233, 1999. 9
Distribution Modeling via Sparse Bayesian Networks
Juan, A. and Vidal, E. On the use of Bernoulli mixture models for text classification. Pattern Recognition, 35 (12):2705–2710, 2002.
Melnykov, V. and Maitra, R. Finite mixture models and model-based clustering. Statistics Surveys, 4:80–116, 2010. doi: 10.1214/09-SS053.
Kneib, T., Silbersdorff, A., and Säfken, B. Rage against the mean–a review of distributional regression approaches. Econometrics and Statistics, 26:99–123, 2023.
Merlevède, F., Peligrad, M., and Rio, E. Bernstein inequality and moderate deviations under strong mixing conditions. In High dimensional probability V: the Luminy volume, volume 5, pp. 273–293. Institute of Mathematical Statistics, 2009.
Kollovieh, M., Ansari, A. F., Bohlke-Schneider, M., Zschiegner, J., Wang, H., and Wang, Y. B. Predict, refine, synthesize: Self-guiding diffusion models for probabilistic time series forecasting. Advances in Neural Information Processing Systems (NeurIPS), 36:28341–28364, 2023.
Morgan, J. N. and Sonquist, J. A. Problems in the analysis of survey data, and a proposal. Journal of the American Statistical Association, 58(302):415–434, 1963. Muzellec, B., Josse, J., Boyer, C., and Cuturi, M. Missing data imputation using optimal transport. In International Conference on Machine Learning (ICML), volume 119, pp. 7130–7140, 2020.
Kotelnikov, A., Baranchuk, D., Rubachev, I., and Babenko, A. Tabddpm: Modelling tabular data with diffusion models. In International Conference on Machine Learning (ICML), pp. 17564–17579, 2023.
Nagler, T. and Czado, C. Evading the curse of dimensionality in nonparametric density estimation with simplified vine copulas. Journal of Multivariate Analysis, 151:69– 89, 2016.
Lee, W. S., Bartlett, P. L., and Williamson, R. C. Efficient agnostic learning of neural networks with bounded fanin. IEEE Transactions on Information Theory, 42(6): 2118–2132, 1996.
Neal, R. M. Connectionist learning of belief networks. Artificial Intelligence, 56(1):71–113, 1992.
Liang, T., Dharmakeerthi, K., and Koriyama, T. Denoising diffusions with optimal transport: Localization, curvature, and multi-scale complexity. Transactions on Machine Learning Research, 2026.
Oh, D. H. and Patton, A. J. Modeling dependence in high dimensions with factor copulas. Journal of Business & Economic Statistics, 35(1):139–154, 2017.
Liu, L., Li, D., and Wong, W. H. Convergence rates of a class of multivariate density estimation methods based on adaptive partitioning. Journal of Machine Learning Research, 24(50):1–64, 2023.
Oko, K., Akiyama, S., and Suzuki, T. Diffusion models are minimax optimal distribution estimators. In International Conference on Machine Learning (ICML), 2023.
Liu, Q. and Paparrizos, J. The elephant in the room: Towards a reliable time-series anomaly detection benchmark. Advances in Neural Information Processing Systems (NeurIPS), 37:108231–108261, 2024.
Prokhorenkova, L., Gusev, G., Vorobev, A., Dorogush, A. V., and Gulin, A. Catboost: Unbiased boosting with categorical features. Advances in Neural Information Processing Systems (NeurIPS), 2018.
Lu, C., Zhou, Y., Bao, F., Chen, J., Li, C., and Zhu, J. Dpmsolver++: Fast solver for guided sampling of diffusion probabilistic models. Machine Intelligence Research, pp. 1–22, 2025.
Rio, E. Covariance inequalities for strongly mixing processes. In Annales de l’IHP Probabilités et statistiques, volume 29, pp. 587–597, 1993. Salimans, T., Karpathy, A., Chen, X., and Kingma, D. P. Pixelcnn++: Improving the pixelcnn with discretized logistic mixture likelihood and other modifications. In International Conference on Learning Representations (ICLR), 2017.
McCracken, M. W. and Ng, S. Fred-md: A monthly database for macroeconomic research. Journal of Business & Economic Statistics, 34(4):574–589, 2016. McElfresh, D., Khandagale, S., Valverde, J., Prasad C, V., Ramakrishnan, G., Goldblum, M., and White, C. When do neural nets outperform boosted trees on tabular data? Advances in Neural Information Processing Systems (NeurIPS), 36:76336–76369, 2023.
Scott, D. W. On optimal and data-based histograms. Biometrika, 66(3):605–610, 1979. Scrucca, L., Fop, M., Murphy, T. B., and Raftery, A. E. mclust 5: clustering, classification and density estimation using Gaussian finite mixture models. The R Journal, 8 (1):289–317, 2016.
Meinshausen, N. and Ridgeway, G. Quantile regression forests. Journal of Machine Learning Research, 7(6), 2006. 10
Distribution Modeling via Sparse Bayesian Networks
Šidák, Z. Rectangular confidence regions for the means of multivariate normal distributions. Journal of the American Statistical Association, 62(318):626–633, 1967. Silveri, M. G. and Ocello, A. Beyond log-concavity and score regularity: Improved convergence bounds for scorebased generative models in w2-distance. In International Conference on Machine Learning (ICML), 2025. Song, Y., Sohl-Dickstein, J., Kingma, D. P., Kumar, A., Ermon, S., and Poole, B. Score-based generative modeling through stochastic differential equations. arXiv preprint arXiv:2011.13456, 2020. Song, Y., JaschaSohl-Dickstein, Kingma, D. P., Kumar, A., Ermon, S., and Poole, B. Score-based generative modeling through stochastic differential equations. In International Conference on Learning Representations (ICLR), 2021. Tauchen, G. Finite state Markov-chain approximations to univariate and vector autoregressions. Economics Letters, 20(2):177–181, 1986. Tibshirani, R. Regression shrinkage and selection via the lasso. Journal of the Royal Statistical Society Series B: Statistical Methodology, 58(1):267–288, 1996. Tsay, R. S. Analysis of financial time series. John wiley & sons, 2005. Tsybakov, A. B. Introduction to Nonparametric Estimation. Springer Series in Statistics. Springer, 2009. Van Den Oord, A., Kalchbrenner, N., and Kavukcuoglu, K. Pixel recurrent neural networks. In International Conference on Machine Learning (ICML), 2016. Vandermeulen, R. A., Tai, W. M., and Aragam, B. Breaking the curse of dimensionality in structured density estimation. Advances in Neural Information Processing Systems, 37:65852–65883, 2024. Zhang, K., Yin, C. H., Liang, F., and Liu, J. Minimax optimality of score-based diffusion models: Beyond the density lower bound assumptions. arXiv preprint arXiv:2402.15602, 2024. Zhang, Z., Hua, Y., Sun, G., Wang, H., and McLoone, S. Training diffusion-based generative models with limited data. In International Conference on Machine Learning (ICML), 2025.
11
Distribution Modeling via Sparse Bayesian Networks
Supplementary Material to “Breaking the Curse with BAND: Nonparametric Distribution Estimation in High Dimensions” This 3.7 and 3.9. Hereafter, we define 3.3 and Corollaries Q Supplementary Material contains the proofs of Theorem #– p p 1{· · · } = 1 when S = ∅, implying P ∩ {· · · } = 1. Let R and R (S, d S ) respectively denote Ac and q∈S c q∈S #– p A(S, d S ) with A = R , as in Section 2.1.
A. Supplementary Material for Main Texts A.1. Empirical Cramér–Wold Distance in Section 4 #– 1 #– 2 Let {X t }nt=1 and {Y t }nt=1 denote two samples of sizes n1 and n2 , respectively. For any unit direction #– v ∈ Rp with #– ∥ v ∥ = 1, define the empirical projected CDFs 2
n1 n #– o 1 X Fb#– 1 #– v⊤X t ≤ z , v (z) = n1 t=1
n2 n #– o 1 X b #– G 1 #– v⊤Y t ≤ z . v (z) = n2 t=1
The empirical Cramér–Wold distance between the two samples is given by b v (z) , Fb#– max max v (z) − G #– #– v ∈V
z∈D
where the set of thresholds D is discretized as D = {−4, −3.9, −3.8, . . . , 7}, to approximate the real line. The set V consists of 10,000 randomly generated unit vectors in Rp , used to approximate the supremum over all directions. #– 1 #– 2 In our experiments, each sequence, {X t }nt=1 or {Y t }nt=1 , corresponds to one of the following: (i) a fitted distribution model, (ii) the ground-truth data-generating process, or (iii) the original training sample. A.2. Monte Carlo Probability Estimation and Sampling via BAND Here, we introduce how to use the fitted BAND model to estimate P(X ∈ A) for any measurable set A ⊆ Rp . Let N denote the number of samples generated from a BAND model. Let #b– #b– { h i,1 , . . . , h i,p }N i=1 #b– denote the generated samples, where each h i,j = (b hi,j,1 , . . . , b hi,j,mj )⊤ is an mj -dimensional one-hot vector. For a fixed sample index i, the sampling procedure proceeds sequentially over coordinates j = 1, . . . , p. For each l ∈ {1, . . . , mj }, define #b– #b– vl = ebjl ( h i,j−1 , . . . , h i,1 ), and apply the thresholding rule #b– #b– eb†jl ( h i,j−1 , . . . , h i,1 ) =
( vl , 0,
mj vl ≥ SPTF × Standard Deviation {vl }l=1 , otherwise.
#b– #b– When j = 1, we set ebjl ( h i,j−1 , . . . , h i,1 ) = ebjl . The ith sample is then assigned to an interval index l ∈ {1, . . . , mj } by sampling from the probability vector #b– #b– #b– #b– mj X eb†jmj ( h i,j−1 , . . . , h i,1 ) eb†j1 ( h i,j−1 , . . . , h i,1 ) #b– #b– ,..., , V = eb†jl ( h i,j−1 , . . . , h i,1 ). V V l=1
#b– Given the sampled index l, we construct the one-hot vector h i,j by setting b hi,j,l ← 1,
b hi,j,k ← 0 12
for all k ̸= l.
Distribution Modeling via Sparse Bayesian Networks
#b– After completing this procedure for all coordinates, we obtain a collection of one-hot vectors { h i,l }1≤i≤N, 1≤l≤p , where #b– #b– each ( h , . . . , h ) identifies a hyper-rectangular cube. For each i, a point #– x is then sampled uniformly from the i,1
i,p
i
corresponding bin. The Monte Carlo estimator of P(X ∈ A) is given by N
1 X #– 1{ x i ∈ A}. N i=1 The above sampling procedure depends on a tuning parameter SPTF ∈ [0, 0.4]. Its optimal value is selected by minimizing the Cramér–Wold distance (see Section 4.2) between the training data and the generated samples. A.3. Conditional Confidence Regions Construction #– The procedure for constructing a k-month forecast confidence region for U t at confidence level α, conditional on the past #– #– state U t−k , using BAND, is as follows. First, BAND selects all hypercubes that contain the past state U t−k ; the predicted probabilities associated with these hypercubes sum to a total probability pb. Second, BAND ranks these hypercubes by their predicted probabilities and selects them sequentially until the cumulative probability reaches αb p. At this point, the union of the selected hypercubes generally overestimates the desired confidence region. To correct this, the final hypercube is shrunk toward its center so that the cumulative probability matches the threshold αb p exactly, under the assumption of uniform probability within each hypercube—a standard assumption for partition-based predictors such as BAND. #– bU The resulting union of hypercubes, denoted by C( t−k ), constitutes the forecasted confidence region. The empirical coverage is then computed as #– #– bU #{t ∈ TestSet | U t ∈ C( t−k )} . #TestSet #– #– When constructing a confidence region for a subset of target variables V t ⊂ U t , the irrelevant coordinates (i.e., those in #– #– U t \ V t ) are marginalized out, and the same procedure is applied. In BAND, this corresponds to summing the predicted probabilities across all hypercubes along the irrelevant dimensions, yielding a confidence region for the variables of interest #– V t. Details of our implementation of BAND are provided in the accompanying R code.
B. Proofs of Main Results B.1. Theorems 3.3 and B.1 and Their Proofs Theorem B.1 below provides the full version of Theorem 3.3, extending it to include results for conditional distributions #– #– while remaining otherwise identical. Let Rpc and Rp (S, d S ) respectively denote Ac and A(S, d S ) with A = Rp . Note that P #– #– #– #– #– the way we construct bins and cubes as in Section 2.2 leads to P(X ∈ A) = #– h 1 ,..., h p P(X ∈ A ∩ Cube( h 1 , . . . , h p )). For the reader’s convenience, we recall the following facts from Sections 2.1–2.2: #– #– #– (1) Cube( h 1 , . . . , h p ) denotes the cube formed by the ωj ( h j )th bin along the jth coordinate; all cubes are disjoint. #– #– #– (2) If #– x ∈ Rp (S, d S ) for some (S, d S ) ∈ Θ, then #– x S lies in the discrete point sets Dj for j ∈ S (i.e., #– x S = d S ), while #– each coordinate of x c is continuous. S
#– #– #– #– #– #– #– (3) For each (S, d S ) ∈ Θ, Cube( h 1 , . . . , h p ) ⊂ Rp (S, d S ) iff Cube( h 1 , . . . , h p ) ∩ Rp (S, d S ) ̸= ∅; additionally, #– #– #– #– Cube( h 1 , . . . , h p ) ⊂ Rpc iff Cube( h 1 , . . . , h p ) ∩ Rpc ̸= ∅. #– #– #– #– (4) Cube( h 1 , . . . , h p ) ⊂ [−δn , δn ]p iff Cube( h 1 , . . . , h p ) ∩ [−δn , δn ]p ̸= ∅. #– #– #– c p (5) πS, #– d S ( x ) = 1 for x in a discrete cube (i.e., S = ∅), and πc ( x ) = 0 when Rc = ∅. #– #– #– Theorem B.1. Let δn > 0 be a positive sequence. √Set ebjl ( h j−1 , . . . , h 1 ) = 0 if ωq ( h q ) = mq for some 1 ≤ q < j or √ T −1 l = mj . For all large n, assume Condition 3.1, ep , and that the pn -dimensional stationary dmin ≤ 1, ε ≤ ( pδn K1 ) 13
Distribution Modeling via Sparse Bayesian Networks
#– distribution of X t satisfies Condition 3.2 with R = δn . It holds for all large n that with probability at least 1 − ∆, for all c #– Q ⊂ {1, . . . , p} with 1 ≤ #Q < p, ϑ > 1, and H ∈ RQ satisfying P(X Qc ∈ H) ≥ ϑ ζ: √ ep T √ #– #– sup P(X ∈ A) − µ b(A) ≤ + ε pδn K1 e + P(X ̸∈ [−δn , δn ]p ) =: ζ, p dmin A∈R and that µ b({ #– x ∈ Rp | #– x Q ∈ B, #– x Qc ∈ H}) 2 #– #– P(X Q ∈ B | X Qc ∈ H) − ≤ . #– p | #– c µ b ({ x ∈ R x (ϑ − 1) ∈ H}) #Q Q B∈R sup
Note that (∆, dmin , T, E † , ε, p) = (∆n , dmin,n , Tn , En† , εn , pn ) but the subscripts n are omitted for simplicity. Proof of Theorem B.1: We begin with defining an auxiliary distribution µ e in (5), which is the population counterpart of #– #– #– our distribution estimator (3). For each ( h 1 , . . . , h p ) ∈ {0, 1}m1 × · · · × {0, 1}mp with one-hot h j , every measurable #– #– A ⊂ Rp , and each #– x ∈ Cube( h 1 , . . . , h p ), define #– #– p X ν̄ #– Y #– #– h 1:p A ∩ Cube( h 1 , . . . , h p ) µ e(A) = ej,ωj (hj ) ( h j−1 , . . . , h 1 ), (5) #– #– ν̄ #– #– #– h 1:p Cube( h 1 , . . . , h p ) j=1 h 1 ,..., h p
#– where ν̄ #– h 1:p denotes Lebesgue measure with respect to the continuous coordinates of Cube( h 1:p ), with zero measure outside #– #– the cube. In addition, if Cube( h 1:p ) is a discrete cube contained by G, then ν̄ #– h 1:p (G) = 1; if Cube( h 1:p ) is a discrete cube not contained by G, then ν̄ #– h 1:p (G) = 0. In what follows, we clarify the definitions of ν̄ #– e, which will be used repeatedly in the sequel. In particular, ν̄ #– h 1:p and µ h 1:p admits a direct integral representation. By fact (3) stated at the beginning of Section B.1, for any A ∈ Rp , ν̄ #– h 1:p (A) satisfies #– #– if A ∩ Cube( h 1 , . . . , h p ) = ∅,
ν̄ #– h 1:p (A) = 0, #– ν̄ #– h 1:p (A) = ν̄ h 1:p
#– A(S, d S ) =
#– ν̄ #– h 1:p (A) = ν̄ h 1:p (Ac ) =
Z #– #– #– #– u ∈A(S, d S )∩Cube( h 1 ,..., h p )
Z #– #– #– u ∈Ac ∩Cube( h 1 ,..., h p )
1 d #– u Sc ,
#– #– #– if A(S, d S ) ∩ Cube( h 1 , . . . , h p ) ̸= ∅,
(6)
#– #– if Ac ∩ Cube( h 1 , . . . , h p ) ̸= ∅.
1 d #– u,
Hereafter, when S c ̸= ∅, integrals with respect to d #– u S c are understood to be taken over the coordinates indexed by S c , #– with the remaining coordinates indexed by S fixed at dRS ; when S c = ∅, such integrals are defined to equal 1. Under #– #– c #– #– 1 d #– this convention, the second integral may be written as #– u S c , where Ā #– h 1:p (S, d S ) = { a S : a ∈ u ∈Ā #h– (S, d S ) 1:p #– #– #– A(S, d S ) ∩ Cube( h 1 , . . . , h p )}. The definition of (6) ensures that all integrals appearing in the subsequent derivations are expressed in a consistent and unified manner. Moreover, by the definition of ν #– h 1:p in (3), #– #– ν̄ #– ν̄ #– h 1:p A ∩ Cube( h 1 , . . . , h p ) h 1:p (A) ν #– #– #– = #– #– . h 1:p (A) = ν̄ #– ν̄ #– h 1:p Cube( h 1 , . . . , h p ) h 1:p Cube( h 1 , . . . , h p ) #– Finally, by the cube construction (following (3)) and the definition of µ e in (5), both P(X ∈ ·) and µ e decompose over the partition: for any A ∈ Rp , X #– #– #– #– P(X ∈ A) = P X ∈ A ∩ Cube( h 1 , . . . , h p ) , #– #– h 1 ,..., h p
µ e(A) =
X #– #– h 1 ,..., h p
#– #– µ e A ∩ Cube( h 1 , . . . , h p ) , 14
(7)
Distribution Modeling via Sparse Bayesian Networks
where the sum ranges over all one-hot vectors. A bias–variance decomposition of the L1 upper bound for our distribution estimator is derived as follows. For any measurable set A ⊂ Rp , #– #– P(X ∈ A) − µ b(A) ≤ µ e(A) − µ b(A) + P(X ∈ A) − µ e(A) . (8) We briefly outline the strategy used to analyze (8). The second term corresponds to the bias. To illustrate how this term is controlled by the estimator defined in (5), we first consider the continuous component while temporarily ignoring tail #– #– #– probabilities. In this case, (5) implies that the density πc ( #– x ) is approximated on each cube by P(X ∈ Cube( h 1 , . . . , h p )), as shown in (10) and (13) below. When the maximum cube diameter ε is sufficiently small, this approximation is accurate provided that πc satisfies the smoothness Condition 3.2. The the variance. Its control relies on Qpfirst term represents #– #– Condition 3.1, which is used to bound the discrepancy between j=1 ej,ωj (hj ) ( h j−1 , . . . , h 1 ) and its empirical counterpart Qp #– #– bj,ωj (hj ) ( h j−1 , . . . , h 1 ), as established in (15) below. j=1 e We now establish an upper bound for the second term on the right-hand side of (8). Analogously to the cube partition argument in (7), we derive a probability upper bound by partitioning according to the discrete–continuous structure of the distribution. In particular, for every measurable set A ⊂ [−δn , δn ]p , we obtain the following partition inequality. #– |P(X ∈ A) − µ e(A)| #– ≤ |P(X ∈ Ac ) − µ e(Ac )| +
X #– (S, d S )∈Θ
#– #– #– |P(X ∈ A(S, d S )) − µ e(A(S, d S ))|.
(9)
We deviation appearing on the right-hand side of (9) individually. To this end, we compare the quantity Qp will bound each #– #– j=1 ej,ωj (hj ) ( h j−1 , . . . , h 1 ) with the corresponding underlying conditional density on each hypercube, as formalized in (12) below. Qp #– #– An expression for j=1 ej,ωj (hj ) ( h j−1 , . . . , h 1 ) is derived as follows. By (1)–(2), together with the definitions of #– the random vectors H j and the probability functions ejl given in Section 2.2, it follows that for every one-hot vector #– h j ∈ {0, 1}mj , p Y
#– #– ej,ωj (hj ) ( h j−1 , . . . , h 1 )
j=1
p \ #– # – = P {H j = h j }
j=1
#– #– #– = P X ∈ Cube( h 1 , . . . , h p ) Z #– #– #– #– c P(X S = d S ) πS, #– dS(u) duS , #– #– #– u ∈Cube( h 1 ,..., h p ) = Z πc ( #– u ) d #– u, #– #– #– u ∈Cube( h 1 ,..., h p )
(10) #– #– #– if Cube( h 1 , . . . , h p ) ∩ Rp (S, d S ) ̸= ∅, #– #– if Cube( h 1 , . . . , h p ) ∩ Rpc ̸= ∅.
Here, the third equality follows from fact (3) stated at the beginning of Section B.1. We also recall the definition of integration with respect to d #– u S c from (6). We next control the discrepancy between these expressions and their approximating densities. By the definition of ε, √ Lemma B.2 in Section B.4, Condition 3.2 with R = δn , and the assumption ε ≤ (K1 pδn )−1 , we obtain that every πθ ⋆ with θ ∈ Θ (see Section 3.1) satisfies sup
∥ #– y ∥2 ≤δn ,∥ #– y − #– x ∥2 ≤ε
√ |πθ ( #– y ) − πθ ( #– x )| ≤ K1 pε × δn × e × πθ ( #– x)
(11)
on ∥ #– x ∥∞ ≤ δn . Combining (11) with the fact that ε denotes the maximum diameter of cubes within [−δn , δn ]p , and 15
Distribution Modeling via Sparse Bayesian Networks
applying the Mean Value Theorem for Integrals, it follows that for each #– x ∈ [−δn , δn ]p , #– #– #– e( #– x )| |P(X S = d S ) πS, #– √ dS(x) − π #– ≤ ε pδn K1 eπS, #– #– #– d S ( x ), P(X S = d S ) √ |πc ( #– x) − π e( #– x )| ≤ ε pδn K1 eπc ( #– x ),
#– #– #– if #– x ∈ Cube( h 1 , . . . , h p ) ∩ Rp (S, d S ), #– #– if #– x ∈ Cube( h 1 , . . . , h p ) ∩ Rpc ,
(12)
which follows from (6) and the definition that for every #– x ∈ [−δn , δn ]p , #– #– R #– #– π #– ( #– u ) d #– u Sc P(X S = d S ) #– u ∈Cube( h 1 ,..., h p ) S, d S , #– #– #– ν̄ Cube( h , . . . , h ) 1 p h 1:p π e( #– x) = R #– #– π ( #– u ) d #– u #– u ∈Cube( h 1 ,..., h p ) c #– , ν̄ #– Cube( #– h 1, . . . , h p) h 1:p
#– #– #– if #– x ∈ Cube( h 1 , . . . , h p ) ∩ Rp (S, d S ), #– #– if #– x ∈ Cube( h 1 , . . . , h p ) ∩ Rpc .
(13)
#– Note that we have defined πc ( #– x ) = 0 for discrete distributions; hence, |πc ( #– x )− π e( #– x )| = 0. Moreover, since πS, #– dS(x) = 1 #– #– #– when Rp (S, d ) consists of a single discrete point, it follows that |P(X = d ) π #– ( #– x) − π e( #– x )| = 0. Therefore, (12) S
S
S
S, d S
may not be tight in this case, but it is sufficient for our analysis.
We now complete the derivation of upper bounds for the terms appearing on the right-hand side of (9). Combining (5)–(7), (10), and (12), together with fact (3) stated at the beginning of Section B.1, we obtain that for every measurable set A ⊂ [−δn , δn ]p , #– |P(X ∈ Ac ) − µ e(Ac )| X #– #– #– #– #– ≤ P(X ∈ Ac ∩ Cube( h 1 , . . . , h p )) − µ e(Ac ∩ Cube( h 1 , . . . , h p )) #– #– h 1 ,..., h p
X
=
#– #– h 1 ,..., h p
≤
Z #– #– #– u ∈Ac ∩Cube( h 1 ,..., h p )
πc ( #– u ) d #– u
Z ν̄ #– h 1:p Ac − πc ( #– u ) d #– u #– #– #– #– #– Cube( h , . . . , h ) ν̄ #– u ∈Cube( h ,..., h ) 1 p 1 p h 1:p X Z |πc ( #– u) − π e( #– u )| d #– u #– #–
#– #– h 1 ,..., h p
#– u ∈Ac ∩Cube( h 1 ,..., h p )
√ ≤ ε pδn K1 e ×
Z #– x ∈Rp c
πc ( #– x )d #– x.
#– Similarly, for each (S, d S ) ∈ Θ, #– #– #– |P(X ∈ A(S, d S )) − µ e(A(S, d S ))| Z X #– #– ≤ P(X S = d S ) #– #– #– h 1 ,..., h p
≤
#– #– #– u ∈A(S, d S ) ∩ Cube( h 1 ,..., h p )
#– #– c πS, #– dS(u) duS
#– Z ν̄ #– #– #– h 1:p A(S, d S ) #– #– c − × P( X = d ) × πS, #– S S #– #– dS(u) duS #– #– #– #– ν̄ h 1:p Cube( h 1 , . . . , h p ) u ∈Cube( h 1 ,..., h p ) X Z #– #– #– P(X S = d S ) × πS, #– e( #– u ) d #– u Sc dS(u) − π #– #– #–
#– #– h 1 ,..., h p
#– u ∈A(S, d S ) ∩ Cube( h 1 ,..., h p )
√ #– #– ≤ ε pδn K1 e × P(X S = d S ) ×
Z #– #– x ∈Rp (S, d S )
16
#– #– c πS, #– d S ( x )d x S .
Distribution Modeling via Sparse Bayesian Networks
With these and (9), we derive that for every measurable A ⊂ [−δn , δn ]p , #– |P(X ∈ A) − µ e(A)| Z √ ≤ ε pδn K1 e × (
#– x ∈Rp c
πc ( #– x )d #– x+
√ ≤ ε pδn K1 e,
X #– (S, d S )∈Θ
#– #– P(X S = d S ) ×
Z #– #– x ∈Rp (S, d S )
#– #– c πS, #– d S ( x )d x S )
(14)
where the third inequality follows from (1). This establishes an upper bound on the second term on the RHS of (8). Now, we turn to establish an upper bound on the first term on the RHS of (8). By Condition 3.1 and the definitions of µ b and µ e respectively in (3) and (5), it holds on E † (given in Condition 3.1), for every measurable A ⊂ [−δn , δn ]p , |e µ(A) − µ b(A)| p Y
X
≤
#– ωq ( h q )<mq for 1≤q≤p j=1
X
Y p
#– ωq ( h q )<mq for 1≤q≤p
j=1
≤
×
p √ X j=1
√ p T ≤ × dmin √ ep T ≤ , dmin
T
dmin
#–
#– #– ebj,ωj ( #– h j ) ( h j−1 , . . . , h 1 )
#– #– ej,ωj ( #– h j ) ( h j−1 , . . . , h 1 ) max
1≤k≤p ωq ( #– z q )<mq for 1≤q≤k
T
p Y j=1
× max √
1+
#–
ej,ωj ( #– h j ) ( h j−1 , . . . , h 1 ) −
p #– #– ebk,ωk ( #– z k ) ( z k−1 , . . . , z 1 ) ∨ 1 #– #– ek,ωk ( #– z k ) ( z k−1 , . . . , z 1 )
(15)
!p
dmin
#– #– #– #– where ωj ( h j ) = q if the qth coordinate of h j equals one, and the summation is taken over all ( h 1 , . . . , h p ) ∈ {0, 1}m1 × #– · · · × {0, 1}mp , with each h j being one-hot and satisfying the specified requirements. Here, we apply the identity for positive ai ’s and bi ’s: ! ! p p j−1 p p Y Y X Y Y ai − bi = (aj − bj ) bi ai i=1
i=1
j=1
=
p Y i=1
i=1
! ai
i=j+1
p X
! j−1 (aj − bj ) Y bi × , aj a j=1 i=1 i
(16)
along with Condition 3.1 and the assumption that #–
#–
min ej,l ( h j−1 , . . . , h 1 ) ≥ dmin , #– j,l<mj ,ωq ( h q )̸=mq for 1≤q<j
(17)
to derive the second inequality. The third inequality follows from Condition 3.1, X
p Y
#– #– h 1 ,..., h p j=1
#– #– ej,ωj (hj ) ( h j−1 , . . . , h 1 ) = 1,
#– #– #– the assumption (17), and the assumption that ebjl ( h j−1 , . . . , h 1 ) = 0 whenever ωq ( h q ) = mq for some 1 ≤ q < j or l = mj . Finally, we use √ √ √ T p T T 1+ = exp p loge 1 + dmin ≤ exp pdmin ≤ e, dmin √
T when pdmin ≤ 1, to obtain the fourth inequality.
17
Distribution Modeling via Sparse Bayesian Networks
By (8), (14)–(15), and the assumption
√
pδn εK1 ≤ 1, we derive that for every measurable A ⊂ [−δn , δn ]p , on E † ,
√ ep T √ #– |P(X ∈ A) − µ b(A)| ≤ ε pδn K1 e + . dmin
(18)
#– #– #– By the definitions of (3) and the assumption that ebjl ( h j−1 , . . . , h 1 ) = 0 whenever ωq ( h q ) = mq for some 1 ≤ q < j or l = mj , we have that µ b(A) = µ b(A ∩ [−δn , δn ]p ). By this, Condition 3.1, and (18), we deduce that, with probability at least 1 − ∆ (on E † ), the following holds for all large n and every A ∈ Rp . #– |P(X ∈ A) − µ b(A)| #– #– ≤ P(X ∈ A\[−δn , δn ]p ) + P(X ∈ A ∩ [−δn , δn ]p ) − µ b(A)
#– #– = P(X ∈ A\[−δn , δn ]p ) + P(X ∈ A ∩ [−δn , δn ]p ) − µ b(A ∩ [−δn , δn ]p ) √ ep T √ #– p ≤ P(X ̸∈ [−δn , δn ] ) + ε pδn K1 e + . dmin
(19)
This completes the proof of the first assertion of Theorem B.1. We now turn to establishing the upper bound for the conditional probability deviation. Recall that
√ ep T √ #– ζ= + ε pδn K1 e + P(X ̸∈ [−δn , δn ]p ). dmin
For any ϑ > 1, we have the following result. Assume a, b > 0 and x, y ∈ (0, 1], with y ≥ ϑζ > ζ, x ≤ y, a ≤ b, and |x − a| ≤ ζ, |y − b| ≤ ζ. Then x a x(b − y) + y(x − a) − = y b yb ≤
|x||b − y| + |y||x − a| |y b|
≤
|x|ζ + |y|ζ |y| |b|
(20)
|x|ζ + |y|ζ ≤ |y| (y − |y − b|) ≤
yζ + yζ y(y − ζ)
≤
2 . ϑ−1
#– #– #– Now, applying the first result of Theorem B.1 and (20) with x = P(X Q ∈ B, X Qc ∈ H), y = P(X Qc ∈ H), a = µ b({ #– x ∈ #– #– #– #– p p b({ x ∈ R | x Qc ∈ H}), we obtain the desired result for the second assertion of R | x Q ∈ B, x Qc ∈ H}), and b = µ Theorem B.1, thereby completing its proof. B.2. Proof of Corollary 3.7 We show that the assumptions of Theorem B.1 are satisfied under the setting of Corollary 3.7. In particular, we verify Condition 3.1 for this setting and then apply Theorem B.1 to conclude the result. We begin by deriving a positive lower bound on dmin > 0 as required by Condition 3.1. 18
Distribution Modeling via Sparse Bayesian Networks
By the definition of ejl in (2) and Condition 3.5(a) or 3.5(b), it holds that #– #– #– #– #– #– ejl ( h j−1 , . . . , h 1 ) = E[H jl | H j−1 = h j−1 , . . . , H 1 = h 1 ] #– #– = E[H jl | H q = h q , q ∈ Sj ] (21) Q #– #– E[H jl × q∈Sj 1{H q = h q }] . = #– #– P(∩q∈Sj {H q = h q }) Q Notation follows Section 2.2, with E[· · · | ∅] = E(· · · ), j∈∅ 1{· · · } = 1, and P(∅) = 1; we will show shortly that the denominator in (21) is strictly positive, so the right-hand side is well defined. In addition, in light of the assumption that the joint density of X J is uniformly bounded above and below by positive constants for every subset J ⊂ {1, . . . , p} with #J ≤ s0 + 1, and the definitions of L1n and L2n , we deduce that there are some positive constants C1 and C2 such that for every 1 ≤ lj < mj , every n ≥ 1, and every J ⊂ {1, . . . , p} with 1 ≤ #J ≤ 1 + s0 , Y (L1n )#J C1 ≤ E H jlj ≤ (L2n )#J C2 . (22) j∈J
By (21)–(22), and the assumptions that #Sj ≤ s0 , L2n ≥ L1n , and supn≥1 L2n L−1 1n is finite, there exists a generic constant C > 0 such that for every n ≥ 1, (L1n )s0 +1 C1 #– #– ejl ( h j−1 , . . . , h 1 ) ≥ (L2n )s0 C2 j,l<mj ,ωq ( h q )̸=mq for 1≤q<j min #–
(23)
≥ L1n C =: dmin .
We now turn to establishing the high probability deviation upper bound as required by the remaining part of Condition 3.1. Analogously to (23), for some generic constant C > 0, it holds for every n ≥ 1 that #–
#–
max ejl ( h j−1 , . . . , h 1 ) ≤ #– j,l<mj ,ωq ( h q )̸=mq for 1≤q<j
(L2n )s0 +1 C2 (L1n )s0 C1
(24)
≤ L2n C. Also, since H jl ∈ {0, 1} almost surely, a direct calculation shows that for every ι ≥ 0, J ⊂ {1, . . . , p}, and 1 ≤ lj ≤ mj , 2+ι Y Y Y E H jlj − E H jlj ≤ E H jlj . (25) j∈J
j∈J
j∈J
#– Here, we recall that H 1jl , . . . , H njl are identically distributed as H jl . We also define sample one-hot H tj = # – (H tj1 , . . . , H tjmj )⊤ , whose stationary population counterpart is denoted by H j = (H j1 , . . . , H jmj )⊤ . Let any #– #– #– j ∈ {1, . . . , p}, l ∈ {1, . . . , mj − 1}, and one-hot ( h 1 , . . . , h p ) with ωq ( h q ) < mq be fixed, and define s Y max1≤q≤p mq #– #– Zt = × H tjl × 1{H tq = h q }. (26) Q #– #– E( q∈Sj 1{H q = h q }) q∈Sj By (22) and (25), it holds that for any 0 < ι ≤ 1, E|Z t − E(Z t )|2+ι ≤ ( max mq )1+ι/2 (L1n )−#Sj − 1≤q≤p
#Sj ×ι 2
−1−ι/2
C1
− 1 (s0 +1)ι
(L2n )#Sj +1 C2 ≤ C5 L1n2
(27)
−1 for some generic constant C5 > 0, since max1≤j≤p mj ≤ L−1 1n C for some constant C > 0 and supn≥1 L2n L1n is finite. −γ1 k By (27), the geometric strong mixing property of {Z t } with α(k) ≤ γ0 e , and the standard covariance inequality for α-mixing sequences (Rio, 1993), we deduce that
2
∞ X
−ι(s0 +1)
|Cov(Z 0 , Z k )| ≤ C × L1n 2+ι
k=0
19
=: v 2
(28)
Distribution Modeling via Sparse Bayesian Networks
for some generic constant C > 0. Here, the mixing property of {Z t } holds because Z t is a measurable function of (X il , l ∈ Sj ∪ {j}) and {(X il , l ∈ Sj ∪ {j})}i is strongly mixing with coefficient α(k) according to Condition 3.6. By Theorem 2 of (Merlevède et al., 2009), (28), that {Z t } is geometrically strongly mixing with α(k) ≤ γ0 e−γ1 k for all k ≥ 1, and that s max1≤j≤p mj −1/2−#Sj /2 |Z t − E(Z t )| ≤ ≤ CL1n =: Qn Q #– #– E( q∈Sj 1{H q = h q }) for some generic constant C > 0 almost surely (derived by using similar calculations for (27)), there exists some C6 > 0 such that for all n ≥ 2, ! n X C 6 x2 . (29) P (Z t − E(Z t )) ≥ x ≤ exp v 2 n + Q2n + xQn (log n)2 t=1 √ −1−#Sj 0 By (29), setting x = v n(log n)4 and using the assumption L1n ≤ L−1−s ≤ n ≤ v 2 n for all large n (the last 1n inequality follows since L1n ≤ 1 for all large n), it follows that, for some constant C > 0 and all sufficiently large n, n
P
−1
n X
! −1/2
(Z t − E(Z t )) ≥ vn
4
(log n)
≤ exp −C(log n)2 .
(30)
t=1
By letting Gt =
q
1 × #– #– P(∩q∈Sj {H q = h q })
Q
#–
#–
q∈Sj 1{H tq = h q } and a similar argument for deriving (30), we obtain some
′
generic constant C > 0 such that for all large n, n
P
−1
n X
! −1/2
(Gt − E(Gt )) ≥ vn
4
(log n)
′ ≤ exp −C (log n)2 .
(31)
t=1
#– #– #– #– Recall the definition of ebjl ( h j−1 , . . . , h 1 ) given in (4). In (33)–(34) below, we establish a bound on |ejl ( h j−1 , . . . , h 1 ) − #– #– ebjl ( h j−1 , . . . , h 1 )|, under the assumption that the complementary events of the left-hand sides of (30) and (31) hold. #– #– #– #– |ejl ( h j−1 , . . . , h 1 ) − ebjl ( h j−1 , . . . , h 1 )| Q Pn Q #– #– #– #– E[H jl × q∈Sj 1{H q = h q }] n−1 t=1 H tjl × q∈Sj 1{H tq = h q } ≤ − Pn Q #– #– #– #– P(∩q∈Sj {H q = h q }) n−1 t=1 q∈Sj 1{H tq = h q } Q Pn Q #– #– #– #– E[H jl × q∈Sj 1{H q = h q }] − n−1 t=1 H tjl × q∈Sj 1{H tq = h q } ≤ #– #– P(∩q∈Sj {H q = h q }) Pn Q #– #– #– #– P(∩q∈Sj {H q = h q }) − n−1 t=1 q∈Sj 1{H tq = h q } #– #– + ebjl ( h j−1 , . . . , h 1 ) × . #– #– P(∩q∈Sj {H q = h q })
(32)
When the complementary event of the left-hand sides of (30) holds, the first term on the RHS of (32) is bounded by s −1/2
vn
4
(log n)
1−#Sj ′ 1 ≤ C vn−1/2 (log n)4 L2n 2 , #– #– P(∩q∈Sj {H q = h q }) × max1≤j≤p mj
′
(33)
for some generic constant C > 0 for all large n, where we use max1≤j≤p mj ≥ L−1 2n C for some constant C > 0, (22), and the assumption that supn≥1 L2n L−1 is finite. In addition, when the complementary events of the left-hand sides of (30) and 1n (31) hold, by the recursive argument together with #– #– #– #– #– #– #– #– ebjl ( h j−1 , . . . , h 1 ) ≤ |ejl ( h j−1 , . . . , h 1 ) − ebjl ( h j−1 , . . . , h 1 )| + ejl ( h j−1 , . . . , h 1 ), 20
Distribution Modeling via Sparse Bayesian Networks
#– #– that 0 ≤ ebjl ( h j−1 , . . . , h 1 ) ≤ 1, (22), (24), and the assumption that supn≥1 L2n L−1 1n is finite, the second term on the RHS of (32) is bounded by s ′ 1 #– #– −1/2 4 C vn (log n) × ejl ( h j−1 , . . . , h 1 ) #– #– P(∩q∈Sj {H q = h q }) s ′′ 1 (34) −1/2 4 ≤ C vn (log n) × L2n #– #– P(∩q∈Sj {H q = h q }) ′′′
1−
#Sj
≤ C vn−1/2 (log n)4 L2n 2 . ′
′′
for positive generic constants C , C , C
′′′
for all large n.
Now, let us deal with the number of deviation bounds needed. The transformed process {Z t } given by (26) is indexed by #– j ∈ {1, . . . , p}, l ∈ {1, . . . , mj − 1}, and ( h q , q ∈ Sj ). As a result, the total number of means to be estimated is bounded by 2 × (m1 +
mj p X X Y
mq ) ≤ 2m1 + 2p( max mj )1+s0 1≤j≤p
j=2 l=1 q∈Sj
≤ 4p( max mj )1+s0
(35)
1≤j≤p
≤ 4n2+s0 , s
− 1 − 0 −b
since p ≤ n (due to our assumption pn−1/2 (log n)5 L2n2 2 0 ≤ 1) and mj ≤ n. The same upper bound holds for the number of {Gt }. Applying (30)–(31) and (35), we conclude that the union of the events corresponding to the left-hand sides of (30)–(31) has negligible probability as n increases; denote this probability by ∆ = ∆n . By this, (28), (30)–(34), and the assumption that L2n = o(1), with √
−ι(s0 +1)
T := CL1n 2+ι
1−s0
n−1/2 (log n)4 L2n2
(36)
for some sufficiently large constant C > 0, we have Condition 3.1 with ∆ = ∆n approaching zero as n increases. In addition, by the assumption that supn≥1 L2n L−1 1n is finite, and (23), it holds for all large n that √ s ep T − 1 − 0 −b ≤ pn−1/2 (log n)5 L2n2 2 0 ≤ 1, dmin
(37)
0 +1) where ι ∈ (0, 1] is arbitrary, s0 is finite, and b0 > 0 is chosen such that b0 = ι×(s . 2+ι
In summary, Condition 3.1 holds for all sufficiently large n, with dmin = dmin,n given by (23) and T = Tn given by (36). √ √ Since ε = εn = L2n p by definition of L2n , the assumption L2n pδn = o(1) implies ε ≤ (δn p K1 )−1 for all large n. Together with (37), the assumption of Condition 3.2 with R = δn , and an application of Theorem B.1, we obtain the desired conclusion of Corollary 3.7. B.3. Proof of Corollary 3.9 We verify that Condition 3.1 holds for all large n in the setting of Corollary 3.9 and then apply Theorem B.1 to establish Corollary 3.7. We begin by deriving a positive constant lower bound dmin > 0 as required by Condition 3.1. By the definition of ejl in (2), the discrete variables assumption, and Condition 3.5(a), it holds that #– #– #– #– #– #– ejl ( h j−1 , . . . , h 1 ) = E[H jl | H j−1 = h j−1 , . . . , H 1 = h 1 ] #– #– = E[H jl | H q = h q , q ∈ Sj ] Q #– #– E[H jl × q∈Sj 1{H q = h q }] . = #– #– P(∩q∈Sj {H q = h q }) 21
(38)
Distribution Modeling via Sparse Bayesian Networks
Q By (38) and the assumption of upper and lower bounds on E( j∈S H jlj ), we deduce that #–
#–
min #– ejl ( h j−1 , . . . , h 1 ) ≥ dmin 1≤j≤p; 1≤l<mj ; ωq ( h q )<mq ,1≤q<j for some constant dmin > 0, which ensures the first part of Condition 3.1. Next, we establish the remaining part of Condition 3.1 for applying Theorem B.1. Let any j ∈ {1, . . . , p}, l ∈ {1, . . . , mj − #– #– #– 1}, and one-hot ( h 1 , . . . , h p ) with ωq ( h q ) < mq be fixed, and define s Y 1 #– #– Zt = 1{H tq = h q }. (39) × H tjl × Q #– #– E( q∈Sj 1{H q = h q }) q∈Sj #– The transformed process {Z t } is indexed by j ∈ {1, . . . , p}, l ∈ {1, . . . , mj − 1}, and ( h q , q ∈ Sj ). Consequently, the total number of distinct transformed processes is bounded as in (35). Q By the assumption of upper and lower bounds on E( j∈S H jlj ) with #S ≤ s0 + 1, max1≤j≤p #Sj ≤ s0 , and that s0 is a constant, it holds that for any 0 < ι ≤ 1, E|Z t − E(Z t )|2+ι ≤ C (40) for some generic constant P C > 0. Moreover, by an argument analogous to that used in (28), we obtain that, for Corollary 3.9, ∞ the upper bound v 2 on 2 k=0 |Cov(Z 0 , Z k )| is constant and that |Z t − E(Z t )| ≤ C =: Qn for some constant C > 0. −1/2 4 Setting x = n (log n) and applying the same arguments as in (29)–(30) and (32)–(34) yields analogous √ √results here (details omitted). Also, the results in (31) carry over directly. Combining these with (35) and T = Tn := Cn−1/2 (log n)4 for a sufficiently large constant C > 0,√we obtain Condition 3.1 with ∆ = ∆n → 0 as n → ∞ for T −1/2 Corollary 3.9 for all large n; additionally, we have that ep (log n)5 ≤ 1 for all large n due to the assumption dmin ≤ pn that pn−1/2 (log n)5 ≤ 1. #– Applying Theorem B.1 with the above result, along with P(X ∈ / [−δn , δn ]) = 0 and ε = 0, yields the desired conclusion of #– #– p p Corollary 3.9. Note that Condition 3.2 is satisfied here because we have defined πS, #– d S ( x ) = 1 on R when R (S, d S ) #– containing a discrete point (see Section 2.1), and that Θ⋆ = Θ when X is discrete (see Section 3.1). We have completed the proof of Corollary 3.9. B.4. Lemma B.2 and its proof √ x such that ∥ #– Lemma B.2. Assume that f ( #– x ) is positive and satisfies ∥∇ log f ( #– x )∥2 ≤ K1 p R for all #– x ∥∞ ≤ R, where #– #– #– #– #– p K > 0 is a constant. Then, for every x , y ∈ R with ∥ x ∥ ≤ R, ∥ y ∥ ≤ R, and ∥ y − #– x ∥ ≤ ε, it holds that ∞
1
∞
√
2
√
|f ( #– y ) − f ( #– x )| ≤ K1 ε pR × f ( #– x ) × exp(K1 pRε), Proof of Lemma B.2: For any continuously differentiable function f : Rp → R and any #– x , #– y ∈ Rp , the multivariate version of the Fundamental Theorem of Calculus (or gradient integral formula along a line) gives Z 1 f ( #– y ) − f ( #– x) = ∇f ( #– x + t( #– y − #– x ))⊤ ( #– y − #– x ) dt. 0
Hence,
|f ( #– y ) − f ( #– x )| ≤ ∥ #– y − #– x ∥2 sup ∥∇f ( #– x + t( #– y − #– x ))∥2 , t∈[0,1]
implying that for every ∥ #– x ∥∞ ≤ R, ∥ #– y ∥∞ ≤ R, ∥ #– x − #– y ∥2 ≤ ε, |f ( #– y ) − f ( #– x )| ≤ ε
sup
∥ #– z ∥∞ ≤R,∥ #– x − #– z ∥2 ≤ε
∥∇f ( #– z )∥2 .
Additionally, we apply the fundamental theorem of calculus along the line segment connecting #– x and #– y: Z 1 ⊤ log f ( #– y ) − log f ( #– x) = ∇ log f #– x + t( #– y − #– x ) ( #– y − #– x ) dt. 0
22
(41)
Distribution Modeling via Sparse Bayesian Networks
√ This, along with our assumption that ∥∇ log f ( #– x )∥2 ≤ K1 p R for all ∥ #– x ∥∞ ≤ R, implies that for any #– x , #– y satisfying #– #– #– #– ∥ x ∥∞ ≤ R, ∥ y ∥∞ ≤ R, and ∥ x − y ∥2 ≤ ε, we have √ | log f ( #– y ) − log f ( #– x )| ≤ ε sup ∥∇ log f ( #– z )∥2 ≤ ε × K1 p R. ∥ #– z ∥∞ ≤R,∥ #– x − #– z ∥2 ≤ε
Hence, for every ∥ #– x ∥∞ ≤ R, ∥ #– y ∥∞ ≤ R, ∥ #– x − #– y ∥2 ≤ ε, f ( #– y) √ exp(log f ( #– y ) − log f ( #– x )) = #– ≤ exp(ε × K1 pR), f(x) leading to that for every ∥ #– x ∥∞ ≤ R, sup
∥ #– z ∥∞ ≤R,∥ #– z − #– x ∥2 ≤ε
√ f ( #– z ) ≤ f ( #– x ) exp(ε × K1 pR).
(42)
By (41)–(42), and our assumption, we deduce that for every ∥ #– x ∥∞ ≤ R, ∥ #– y ∥∞ ≤ R, ∥ #– x − #– y ∥2 ≤ ε, |f ( #– y ) − f ( #– x )| ≤ ε =ε
sup
∥∇f ( #– z )∥2
sup
[∥∇ log f ( #– z )∥2 f ( #– z )]
∥ #– z ∥∞ ≤R,∥ #– z − #– x ∥2 ≤ε ∥ #– z ∥∞ ≤R,∥ #– z − #– x ∥2 ≤ε
≤ ε × [ sup
∥ #– z ∥∞ ≤R
∥∇ log f ( #– z )∥2 ] ×
sup
∥ #– z ∥∞ ≤R,∥ #– z − #– x ∥2 ≤ε
f ( #– z)
√ √ ≤ ε × K1 pR × f ( #– x ) × exp(K1 pRε),
where the first inequality follows from (41), the first equality is due to f ( #– x ) > 0 on ∥ #– x ∥∞ ≤ R, and the third inequality results from our assumptions.
C. Supplementary Material D.1. Examples of Distributions Satisfying Condition 3.2 We present three examples of distributions satisfying Condition 3.2. To handle the non-differentiability of the uniform distribution at its boundary, we consider a slightly enlarged support in Example C.3 (or, alternatively, a smoothly clipped version). Note that finite mixtures of these distributions, such as Gaussian or Student-t mixtures, similarly satisfy Condition 3.2, though we omit those derivations for brevity. Example C.1. Let f ( #– x ) be the multivariate centered Student-t density with a uniformly positive definite scale matrix Σ ∈ Rp×p and degrees of freedom ν > 0. Assume p/ν is uniformly bounded. Then, Condition 3.2 holds for every R > 0. Example C.2. Let f ( #– x ) be the multivariate Gaussian distribution with a uniformly positive definite covariance matrix. Then, Condition 3.2 holds for every R > 0. Example C.3. Let f ( #– x ) be a uniform distribution on [−δn − ϵ, δn + ϵ]p for some small ϵ > 0. Then, Condition 3.2 holds for R = δn . In Examples C.1–C.3, constant K1 of Condition 3.2 is independent of R. Proof of Examples C.1–C.3: We begin with establishing Condition 3.2 for Example C.1. Let f be the multivariate Student-t density with location #– µ ∈ Rp , positive definite scale matrix Σ ∈ Rp×p , and degrees of freedom ν > 0. Specifically, −(ν+p)/2 Γ ν+p 1 #– #– ⊤ −1 #– #– 2 f ( #– x) = ν 1 + ( x − µ ) Σ ( x − µ ) , ν Γ 2 (νπ)p/2 D(Σ)1/2 µ is where D(Σ) denotes the determinant of the scale (or covariance) matrix. Additionally, we assume in Example C.1 that #– a vector of zeros. Let us begin the formal proof. The derivative of log-density is given by ν+p Σ−1 #– x ∇ log f ( #– x) = − × . #– ⊤ −1 ν 1 + x Σ #– x /ν 23
Distribution Modeling via Sparse Bayesian Networks
Hence
x ∥2 ν+p ν + p ∥Σ−1 ∥2 ∥ #– ≤ ∥Σ−1 ∥2 ∥ #– x ∥2 , ∥∇ log f ( #– x )∥2 ≤ ν 1 + #– x ⊤ Σ−1 #– x /ν ν √ since (1 + #– x ⊤ Σ−1 #– x /ν)−1 ≤ 1. If max |x | ≤ R then ∥ #– x ∥ ≤ p R, so j
j
2
ν+p √ ∥∇ log f ( #– x )∥2 ≤ ∥Σ−1 ∥2 p R. ν
(D.1)
Thus one can take K1 = (ν + p)/ν ∥Σ−1 ∥2 to conclude Condition 3.2 for every R > 0. We have completed the proof of Example C.1. Now, we turn to establish Condition 3.2 for Example C.2. Let f denote the multivariate Gaussian density: 1 ⊤ −1 #– #– #– 1 #– exp − f ( #– x) = ( x − µ ) Σ ( x − µ ) , 2 (2π)p/2 D(Σ)1/2 where #– µ ∈ Rp is the mean vector, and Σ ∈ Rp×p is a positive definite covariance matrix. In Example C.2, we consider the scenario with #– µ being a vector of zeros. The log-density is given by p 1 1 ⊤ −1 #– log f ( #– x ) = − log(2π) − log D(Σ) − #– x Σ x. 2 2 2 Hence, the gradient of the log-density is ∇ log f ( #– x ) = −Σ−1 #– x. Thus, for #– x ∈ Rp with ∥ #– x ∥∞ ≤ R, √ ∥∇ log f ( #– x )∥2 = ∥Σ−1 #– x ∥2 ≤ ∥Σ−1 ∥2 ∥ #– x ∥2 ≤ ∥Σ−1 ∥2 pR. Using this result and setting K1 = ∥Σ−1 ∥2 , we establish Condition 3.2 for all R > 0 in Example C.2. Meanwhile, the analysis for Example C.3 is straightforward and therefore omitted. This completes the proofs for Examples C.1–C.3. D.2. Proof for Section 4.1 When r = 0, Condition 3.5(a) is satisfied under Examples 4.1–4.2 because X 1 , . . . , X p are i.i.d. On the other hand, consider Example 4.2 with r > 3. If we set δn ≥ 32 + r and choose sufficiently many evenly spaced splitting points, then for each variable j ∈ {1, . . . , p}, every one-dimensional bin interval B ⊂ R satisfies either sup B ≤ 1.5 or inf B ≥ 1.5 (bins are mutually disjoint). Hence each bin intersects with at most one of the two modes, which implies that #– the mode of the jth variable is already determined once H l is observed for any 1 ≤ l < j in this scenario. Therefore, for every j > 1, each 1 ≤ k ≤ mj , and each l ∈ {j − 1, . . . , 1}, #– #– #– P H jk = 1 H j−1 , . . . , H 1 = P H jk = 1 H l , and Condition 3.5(a) holds in this case with Sj = {l} for any 1 ≤ l < j. D.3. Strong Mixing Processes #– A strong mixing or α-mixing process is defined as follows. A stochastic process {X t }t≥1 is said to be strongly mixing (or α-mixing) if there exists a sequence {α(k)}k≥1 with α(k) → 0 as k → ∞ such that α(k) := sup
sup
∞ t≥1 A∈F1t , B∈Ft+k
P(A ∩ B) − P(A)P(B) ,
#– #– ∞ where F1t and Ft+k are the σ-fields generated by {X s : 1 ≤ s ≤ t} and {X s : s ≥ t + k}, respectively. If there exist constants γ0 , γ1 > 0 such that α(k) ≤ γ0 e−γ1 k ∀k ≥ 1, the process is called geometrically strongly mixing. 24