--- title: "ROME: Robust Multi-Modal Density Estimator" authors: - "Anna Mészáros" - "Julian F. Schumann" - "Javier Alonso-Mora" - "Arkady Zgonnikov" - "Jens Kober" doi: "10.24963/ijcai.2024/525" url: "https://ijcai.org/proceedings/2024/525" pdf: "https://ijcai.org/proceedings/2024/0525.pdf" venue: "Proceedings of the Thirty-Third International Joint Conference on Artificial Intelligence" year: 2024 pages: "4751-4759" track: "Main Track" keywords: - "Machine Learning: ML: Probabilistic machine learning" - "Machine Learning: ML: Evaluation" --- # ROME: Robust Multi-Modal Density Estimator Anna Mészáros, Julian F. Schumann, Javier Alonso-Mora, Arkady Zgonnikov, Jens Kober *Proceedings of the Thirty-Third International Joint Conference on Artificial Intelligence*, pages 4751-4759. DOI: [10.24963/ijcai.2024/525](https://doi.org/10.24963/ijcai.2024/525) ## Abstract The estimation of probability density functions is a fundamental problem in science and engineering. However, common methods such as kernel density estimation (KDE) have been demonstrated to lack robustness, while more complex methods have not been evaluated in multi-modal estimation problems. In this paper, we present ROME (RObust Multi-modal Estimator), a non-parametric approach for density estimation which addresses the challenge of estimating multi-modal, non-normal, and highly correlated distributions. ROME utilizes clustering to segment a multi-modal set of samples into multiple uni-modal ones and then combines simple KDE estimates obtained for individual clusters in a single multi-modal estimate. We compared our approach to state-of-the-art methods for density estimation as well as ablations of ROME, showing that it not only outperforms established methods but is also more robust to a variety of distributions. Our results demonstrate that ROME can overcome the issues of over-fitting and over-smoothing exhibited by other estimators. ## Full text ## 1 Introduction Numerous processes are non-deterministic by nature, from geological and meteorological occurrences, biological activities, as well as the behavior of living beings. Estimating the underlying probability density functions (PDFs) of such processes enables a better understanding of them and opens possibilities for probabilistic inference regarding future developments. Density estimation is instrumental in many applications including classification, anomaly detection, and evaluating probabilistic AI models, such as generative adversarial networks [Goodfellow et al., 2020], variational autoencoders [Wei et al., 2020], normalizing flows [Kobyzev et al., 2020], and their numerous variations. When probabilistic models are trained on multi-modal data, they are often evaluated using simplistic metrics (e.g., mean squared error (MSE) between the predicted and ground truth samples). However, such simplistic metrics are unsuited for determining how well a predicted distribution corresponds to the underlying distribution, as they do not capture the fit to the whole distribution. For example, the lowest MSE value between true and predicted samples could be achieved by accurate predictions of the mean of the true underlying distribution whereas potential differences in variance or shape of the distribution would not be penalized. This necessitates more advanced metrics that evaluate the match between the model and the (potentially multi-modal) data. For instance, negative log-likelihood (NLL), Jensen-Shannon divergence (JSD), and expected calibration error (ECE) can be used to evaluate how well the full distribution of the data is captured by the learned models [Xu et al., 2019; Mozaffari et al., 2020; Rasouli, 2020]. However, most data-driven models represent the learned distribution implicitly, only providing individual samples and not the full distribution as an output. This complicates the comparison of the model output to the groundtruth data distributions since the above metrics require distributions, not samples, as an input. Practically, this can be addressed by estimating the predicted probability density based on samples generated by the model. ∗Extended version with appendix: https://arxiv.org/abs/2401.10566 Figure 1: ROME takes samples from unknown distributions and estimates their densities to enable further downstream applications. ![Image](https://ijcai.org/proceedings/2024/fulltext/525/image_000000_bde47270c2afe0702508defe85b4a9949251cd94e4de9876d0e394a9aaa02006.png) Simple methods like Gaussian mixture models (GMM), kernel density estimation (KDE) [Deisenroth et al., 2020], and k-nearest neighbors (kNN) [Loftsgaarden and Quesenberry, 1965] are commonly used for estimating probability density functions. These estimators however rely on strong assumptions about the underlying distribution, and can thereby introduce biases or inefficiencies in the estimation process. For example, problems can arise when encountering multi-modal, non-normal, and highly correlated distributions (see Section 2). While more advanced methods such as manifold Parzen windows (MPW) [Vincent and Bengio, 2002] and vine copulas (VC) [Nagler and Czado, 2016] exist, they have not been thoroughly tested on such problems, which raises questions about their performance. To overcome these limitations, we propose a novel density estimation approach: RObust Multi-modal Estimator (ROME). ROME employs a non-parametric clustering approach to segment potentially multi-modal distributions into separate uni-modal ones (Figure 1). These uni-modal subdistributions are then estimated using a downstream probability density estimator (such as KDE). We test our proposed approach against a number of existing density estimators in three simple two-dimensional benchmarks designed to evaluate a model's ability to successfully reproduce multi-modal, highly-correlated, and non-normal distributions. Finally, we test our approach in a real-world setting using a distribution over future trajectories of human pedestrians created based on the Forking Paths dataset [Liang et al., 2020]. ## 2 RelatedWork The most common class of density estimators are so-called Parzen windows [Parzen, 1962], which estimate the density through an aggregation of parametric probability density functions. A number of common methods use this approach, with KDE being a common non-parametric method [Silverman, 1998]. Provided a type of kernel - which is oftentimes a Gaussian but can be any other type of distribution - KDE places the kernels around each sample of a data distribution and then sums over these kernels to get the final density estimation over the data. This method is often chosen as it does not assume any underlying distribution type [Silverman, 1998]. However, if the underlying distribution is highly correlated, then the common use of a single isotropic kernel function can lead to over-smoothing [Vincent and Bengio, 2002; Wang et al., 2009]. Among multiple approaches for overcoming this issue [Wang et al., 2009; Gao et al., 2022], especially noteworthy is the MPW approach [Vincent and Bengio, 2002]. It uses a unique anisotropic kernel for every datapoint, estimated based on the correlation of the k-nearest neighbors of each sample. However, it has not been previously tested in high-dimensional benchmarks, which is especially problematic as the required memory scales quadratically with the dimensionality of the problem. Another common subtype of Parzen windows are GMMs [Deisenroth et al., 2020], which assume that the data distribution can be captured through a weighted sum of Gaussian distributions. The parameters of the Gaussian distributions - also referred to as components - are estimated through expected likelihood maximization. Nonetheless, especially for non-normal distributions, one needs to have prior knowledge of the expected number of components to achieve a good fit without over-fitting to the training data or oversmoothing [McLachlan and Rathnayake, 2014]. Besides different types of Parzen windows, a prominent alternative is kNN [Loftsgaarden and Quesenberry, 1965] which uses the local density of the k nearest neighbors of every data point to estimate the overall density. While this method is non-parametric, it cannot be guaranteed that the resulting distribution will be normalized [Silverman, 1998]. This could be rectified by using Monte Carlo sampling to obtain a good coverage of the function's domain and obtain an accurate estimate of the normalization factor, which, however, becomes computationally intractable for high-dimensional distributions. When it comes to estimating densities characterized by correlations between dimensions, copula-based methods are an often favored approach. Copula-based methods decompose a distribution into its marginals and an additional function, called a copula, which describes the dependence between these marginals over a marginally uniformly distributed space. The downside of most copula-based approaches is that they rely on choosing an adequate copula function (e.g., Gaussian, Student, or Gumbel) and estimating their respective parameters [Joe, 2014]. One non-parametric copula-based density estimator [Bakam and Pommeret, 2023] aims to address this limitation by estimating copulas with the help of superimposed Legendre polynomials. While this can achieve good results in estimating the density function, it may become computationally intractable as the distribution's dimensionality increases. Another approach involves the use of VC [Nagler and Czado, 2016], which assume that the whole distribution can be described as a product of bivariate Gaussian distributions, thus alleviating the curse of dimensionality. Its convergence, however, can only be guaranteed for normal distributions. Elsewhere [Otneim and Tjøstheim, 2017], a similar approach was pursued, with changes such as using logarithmic splines instead of univariate KDEs for estimating the marginals. However, both of these approaches are not designed for multi-modal distributions and have not been thoroughly tested on such problems. ## 3 RObust Multi-modal Estimator (ROME) The problem of density estimation can be formalized as finding a queryable ̂ p ∈ P, where such that ̂ p is close to the non-queryable PDF p underlying the N available M-dimensional samples X ∈ RN×M : X ∼ p. A solution to the above problem would be an estimator f : RN×M →P, resulting in ̂ p = f(X). Our proposed estimator fROME 1 (Algorithm 1) is built on top of non-parametric cluster extraction. Namely, by separating groups of samples surrounded by areas of low density - expressing the mode of the underlying distribution - we reduce the multi-modal [1 Source code: https://github.com/anna-meszaros/ROME](https://github.com/anna-meszaros/ROME) ## Algorithm 1 ROME ## function TRAINROME(X) density estimation problem to multiple uni-modal density estimation problems for each cluster. The distributions within each cluster then become less varied in density or correlation than the full distribution. Combining this with decorrelation and normalization, the use of established methods such as KDE to estimate probability densities for those unimodal distributions is now more promising, as problems with multi-modality and correlated modes (see Section 2) are accounted for. The multi-modal distribution is then obtained as a weighted average of the estimated uni-modal distributions. ## 3.1 Extracting Clusters To cluster samples X, ROME employs the OPTICS algorithm [Ankerst et al., 1999] that can detect clusters of any shape with varying density using a combination of reachability analysis - which orders the data in accordance to reachability distances - followed by clustering the ordered data based on these reachability distances. In the first part of the algorithm, the reachability analysis is used to sequentially transfer samples from a set of unincluded samples XU,i to the set of included and ordered samples XI,i, starting with a random sample x1 (XI,1 = {x1} and XU,1 = X\{x1}). The samples xi+1 are then selected at iteration i based on the reachability distance r: This sample is then transferred between sets, with XI,i+1 = XI,i ∪ {xi+1} and XU,i+1 = XU,i\{xi+1}, while expanding the reachability set Ri+1 = Ri ∪ {r(xi+1,XI,i)} (with R1 = {∞}). The reachability distance dr in Equation (1) is defined as where minkc is the kc-smallest value of all the available ̂x, used to smooth out random local density fluctuations. We use where kmin and kmax ensure there are sufficient but not too many points for this smoothing, while the term NM/αk adjusts kc to the number of samples and dimensions. The second part of the OPTICS algorithm - after obtaining the reachability distances RN and the ordered set XI,N - is the extraction of a set of clusters C, with clusters C = {cmin, . . . , cmax} ∈ C (with XC = {xI,N,j | j ∈ C} ∈ R|C|×M ). As the computational cost of creating such a cluster set is negligible compared to the reachability analysis, we can test multiple clusterings generated using two different approaches (with rbound = min{rN,cmin , rN,cmax+1}, see Appendix A for further discussion): - First, we use DBSCAN [Ester et al., 1996] for generating the clustering Cϵ based on an absolute limit ϵ, where a cluster must fulfill the condition: - Second, we use ξ-clustering [Ankerst et al., 1999] to generate the clustering Cξ based on the proportional limit ξ, where a cluster fulfills: In both cases, each prospective cluster C also has to fulfill the condition |C| ≥ 2. However, it is possible that not every sample can be fit into a cluster fulfilling the conditions above. These samples are then kept in a separate noise cluster Cnoise that does not have to fulfill those conditions (Cnoise ∈ Cϵ or Cnoise ∈ Cξ respectively). Upon generating multiple sets of clusters Cϵ (ϵ ∈ ε) and Cξ (ξ ∈ ξ), we select the final set of clusters C that achieves the highest silhouette score 2 S = SIL(C,XI,N) ∈ [-1, 1] [Rousseeuw, 1987]. The clustering then allows us to use PDF estimation methods on uni-modal distributions. 2 The silhouette score measures the similarity of each object to its own cluster's objects compared to the other clusters' objects. Figure 2: Samples from the two-dimensional synthetic distributions used for evaluating different PDF estimators. ![Image](https://ijcai.org/proceedings/2024/fulltext/525/image_000001_0aca5c7d5220349be64b4dc8b887112eff54335d6be74053996d147735229101.png) ## 3.2 Feature Decorrelation In much of real-life data, such as the distributions of a person's movement trajectories, certain dimensions of the data are likely to be highly correlated. Therefore, the features in each cluster C ∈ C should be decorrelated using a rotation matrix RC ∈ RM×M . In ROME, RC is found using principal component analysis (PCA) [Wold et al., 1987] on the cluster's centered samples XC = XC - xC (xC is the cluster's mean value). An exception are the noise samples in Cnoise, which are not decorrelated (i.e., RCnoise = I). One can then get the decorrelated samples XPCA,C: ## 3.3 Normalization After decorrelation, we use the matrix ˜ΣC ∈ RM×M to normalize XPCA,C : Here, ˜ΣC is a diagonal matrix with the entries ˜ σ(m) C . To avoid over-fitting to highly correlated distributions, we introduce a regularization with a value σmin (similar to [Vincent and Bengio, 2002]) that is applied to the empirical standard deviation σ(m) PCA,C = √Vm (XPCA,C) (Vm is the variance along feature m) of each rotated feature m (with C′ = C\{Cnoise}): . ## 3.4 Estimating the Probability Density Function Taking the transformed data (decorrelated and normalized), ROME fits the Gaussian KDE fKDE on each separate cluster C as well as the noise samples Cnoise. For a given bandwidth bC for data in cluster C, this results in a partial PDF ̂ pC . The bandwidth bC is set using Silverman's rule [Silverman, 1998]: Figure 3: Samples from the multi-modal pedestrian trajectory distribution [Liang et al., 2020] used for evaluating different PDF estimators. Trajectories span 12 timesteps recorded at 2.5Hz. ![Image](https://ijcai.org/proceedings/2024/fulltext/525/image_000002_e2b41705a0182f69cbe51743990a8d46648726be166445cc84cd7e2dbbafa84a.png) To evaluate the density function ̂ p = fROME(X) for a given sample x, we take the weighted averages of each cluster's ̂ pC : Here, the term |C|/N is used to weigh the different distributions of each cluster with the size of each cluster, so that each sample is represented equally. As the different KDEs ̂ pC are fitted to the transformed samples, we apply them not to the original sample x, but instead apply the identical transformation used to generate those transformed samples, using ̂ pC ((x-xC)). To account for the change in density within a cluster C introduced by the transformation TC, we use the factor | det (TC) |. ## 4 Experiments We compare our approach against two baselines from the literature (VC [Nagler and Czado, 2016] and MPW [Vincent and Bengio, 2002]) in four scenarios, using three metrics. Additionally, we carry out an ablation study on our proposed method ROME. For the hyperparameters pertaining to the clustering within ROME (see Section 3.1), we found empirically that stable results can be obtained using 199 possible clusterings, 100 for DBSCAN (Equation (3)) combined with 99 for ξ-clustering (Equation (4)) as well as using kmin = 5, kmax = 20, and αk = 400 for calculating kc (Equation (2) and Appendix C). ## 4.1 Distributions In order to evaluate different aspects of a density estimation method f , we used a number of different distributions. Table 1: Baseline Comparison - marked in red are cases with notably poor performance; best values are underlined. | Distrib. | ROME | DJS ↓10 MPW | VC | ROME | ̂ W → 0 MPW | VC | ROME | ̂L ↑ MPW | VC | |--------------|-------------|---------------|-------------|------------|---------------|-----------|------------|------------|--------------| | Aniso | 0.010±0.001 | 0.026±0.002 | 0.005±0.001 | -0.13±0.31 | -0.60±0.13 | 1.91±0.91 | -2.53±0.02 | -2.57±0.02 | -3.19±0.02 | | Varied | 0.011±0.001 | 0.025±0.002 | 0.008±0.001 | -0.13±0.20 | -0.49±0.11 | 1.27±0.53 | -4.10±0.03 | -4.12±0.03 | -4.29±0.03 | | Two Moons | 0.002±0.001 | 0.023±0.002 | 0.008±0.002 | 1.40±0.52 | -0.52±0.10 | 1.36±0.51 | -1.02±0.01 | -0.36±0.01 | -0.95±0.01 | | Trajectories | 0.008±0.002 | 0.016±0.001 | 0.743±0.005 | 1.03±0.22 | 1.24±0.24 | 9.30±1.30 | 29.32±0.02 | 26.09±0.02 | -215.23±17.6 | Figure 4: Samples obtained with ROME, MPW and VC (pink) contrasted with samples from p (blue); Varied. ![Image](https://ijcai.org/proceedings/2024/fulltext/525/image_000003_b71096e7c79ded4ce5c37ace09d953e8a6775156bb21321ebed206215c4f0fb3.png) - Three two-dimensional synthetic distributions (Figure 2) were used to test the estimation of distributions with multiple clusters, which might be highly correlated (Aniso) or of varying densities (Varied), or express nonnormal distributions (Two Moons). - A multivariate, 24-dimensional, and highly correlated distribution generated from a subset of the Forking Paths dataset [Liang et al., 2020] (Figure 3). The 24 dimensions correspond to the x and y positions of a human pedestrian across 12 timesteps. Based on 6 original trajectories (x∗i ∈ R12×2 ), we defined the underlying distribution p in such a way, that one could calculate a sample x ∼ p with: Here, Rθ ∈ R2×2 is a rotation matrix rotating x∗i by θ ∼N(0, π 180 ), while s ∼ N (1, 0.03) is a scaling factor. n = N(0, 0.03I) ∈ R12×2 is additional noise added on all dimensions using L ∈ R12×12 , a lower triangular matrix that only contains ones. Further tests on uni-modal problems can be found in Appendix E. ## 4.2 Evaluation andMetrics When estimating density ̂ p, since we cannot query the distribution p underlying the samples X, we require metrics that can provide insights purely based on those samples. To this end we use the following three metrics to quantify how well a given density estimator f can avoid both over-fitting and over-smoothing. - To test for over-fitting, we first sample two different datasets X1 and X2 (N samples each) from p (X1,X2 ∼ p). We then use the estimator f to create two queryable distributions ̂ p1 = f(X1) and ̂ p2 = - f(X2). If those distributions ̂ p1 and ̂ p2 are similar, it would mean the tested estimator does not over-fit; we measure this similarity using the Jensen-Shannon divergence [Lin, 1991]: This metric, however, is not able to account for systematic biases that could be present in the estimator f . Comparisons of ̂ p1 with p - if the dataset allows - can be found in Appendix D. - To test the goodness-of-fit of the estimated density, we first generate a third set of samples ̂ X1 ∼ ̂ p1 with N samples. We then use the Wasserstein distance W [Villani, 2009] on the data to calculate the indicator ̂ W: Here, ̂ W > 0 indicates over-smoothing or misplaced modes, while -1 ≤ ̂ W < 0 indicates over-fitting. - Not every density estimator f has the ability to generate the samples ̂ X1. Consequently, we need to test for goodness-of-fit without relying on ̂ X1. Therefore, we use the average log-likelihood which would be maximized only for p = ̂ p1 as long as X2 is truly representative of p (see Gibbs' inequality [Kvalseth, 1997]). However, using this metric might be meaningless if ̂ p1 is not normalized, as the presence of the unknown scaling factor makes the ̂L values of two estimators incomparable, and cannot discriminate between over-fitting and over-smoothing. For each candidate method f, we used N = 3000, and every metric calculation was repeated 100 times to take into account the inherent randomness of sampling from the distributions, with the standard deviation in the tables being reported with ±. Table 2: Ablations (̂ W → 0, Varied): Clustering is essential to prevent over-smoothing. ROME highlighted in gray. Note that the differences between Silhouette and DBCV are not statistically significant. | Cluster. | Norm. - Decorr. | Norm. - No decorr. | No norm. - No decorr. | fGMM | |------------|-------------------|----------------------|-------------------------|------------| | Silhouette | -0.13±0.20 | -0.13±0.20 | -0.13±0.19 | -0.14±0.20 | | DBCV | -0.09±0.23 | -0.09±0.23 | -0.07±0.23 | -0.03±0.24 | | No clus. | 2.28±0.72 | 2.26±0.72 | -0.17±0.26 | 10.53±2.54 | ## 4.3 Ablations To better understand the performance of our approach, we investigated variations in four key aspects of ROME: - Clustering approach. First, we replaced the silhouette score (see Section 3.1) with density based cluster validation (DBCV) [Moulavi et al., 2014] when selecting the optimal clustering out of the 199 possibilities. Furthermore, we investigated the approach of no clustering (C = {{1, . . . ,N}}) . - Decorrelation vs No decorrelation. We investigated the effect of removing rotation by setting RC = I. - Normalization vs No normalization. We studied the sensitivity of our approach to normalization by setting ˜ΣC = I. - Downstream density estimator. We replaced fKDE with two other candidate methods. First, we used a singlecomponent Gaussian mixture model fGMM fitted to the observed mean ̂µ and covariance matrix ̂Σ of a dataset X. Second, we used a k-nearest neighbor approach fkNN [Loftsgaarden and Quesenberry, 1965] where VM is the volume of the M-dimensional unit hypersphere. We used the rule-of-thumb k = ⌊√N⌋. However, this estimator cannot generate samples. While those four factors would theoretically lead to 24 estimators, fKDE as well as fkNN being invariant against rotation and fGMM being invariant against any linear transformation means that only 14 of ROME's ablations are actually unique. ## 5 Results ## 5.1 Baseline Comparison We found that ROME avoids the major pitfalls displayed by the two baseline methods on the four tested distributions (Table 1). Out of the two baseline methods, the manifold Parzen windows (MPW) approach has a stronger tendency to overfit in the case of the two-dimensional distributions compared to ROME, as quantified by lower DJS values achieved by ROME. MPW does, however, achieve a better log-likelihood for the Two Moons distribution compared to ROME. This could be due to the locally adaptive non-Gaussian distributions being less susceptible to over-smoothing than our approach of using a single isotropic kernel for each cluster if such clusters are highly non-normal. Lastly, in the case of the pedestrian trajectories distribution, ROME once more achieves better performance than MPW, with MPW performing worse both in terms of DJS and ̂L. Table 3: Ablations (̂L ↑, Aniso): When clustering, decorrelation and normalization improve results for distributions with high intra-mode correlation. ROME highlighted in gray. | Cluster. | Norm. Decorr. | No decorr. | No norm. | |------------|-----------------|--------------|------------| | Silhouette | -2.53±0.02 | -2.70±0.01 | -2.79±0.01 | | DBCV | -2.56±0.02 | -2.69±0.01 | -2.83±0.02 | Table 4: Ablations (̂ W → 0, Two Moons): Excluding normalization or using fGMM as the downstream estimator is not robust against nonnormal distributions. ROME highlighted in gray. | Cluster. | Norm. - Decorr. | Norm. - No decorr. | No norm. - No decorr. | fGMM | |------------|-------------------|----------------------|-------------------------|-----------| | Silhouette | 1.40±0.52 | 1.41±0.52 | 3.65±0.99 | 4.37±1.14 | | DBCV | 1.59±0.56 | 1.59±0.56 | 4.24±1.13 | 4.77±1.23 | | No clus. | 1.82±0.60 | 1.84±0.60 | 3.17±0.89 | 5.29±1.34 | Meanwhile, the vine copulas (VC) approach tends to oversmooth the estimated densities (large positive ̂ W values), and even struggles with capturing the different modes (see Figure 4). This is likely because VC uses KDE with Silverman's rule of thumb, which is known to lead to over-smoothing in the case of multi-modal distributions [Heidenreich et al., 2013]. Furthermore, on the pedestrian trajectories distribution, we observed both high DJS and ̂ W values, indicating that VC is unable to estimate the underlying density; this is also indicated by the poor log-likelihood estimates. Overall, while the baselines were able to achieve better performance in selected cases (e.g., MPW better than ROME in terms of ̂ W and ̂L on the Two Moons distribution), they have their apparent weaknesses. Specifically, MPW achieves poor results for most metrics in the case of varying densities within the modes (Varied), while Vine Copulas obtain the worst performance across all three metrics in the case of the multivariate trajectory distributions. ROME, in contrast, achieved high performance across all the test cases. ## 5.2 Ablation Studies When it comes to the choice of the clustering method, our experiments show no clear advantage for using either the silhouette score or DBCV. But as the silhouette score is computationally more efficient than DBCV, it is the preferred method. Table 5: Ablations (DJS ↓10, Trajectories; values are multiplied by 10 for easier comprehension): Using fkNN as the downstream estimator tends to lead to over-fitting. ROME highlighted in gray. | Clustering | Normalization - Decorrelation - fKDE | Normalization - Decorrelation - fkNN | Normalization - No decorrelation - fKDE | No decorrelation - fkNN | No normalization - No decorrelation - fKDE | No normalization - No decorrelation - fkNN | fGMM | |--------------|----------------------------------------|----------------------------------------|-------------------------------------------|---------------------------|----------------------------------------------|----------------------------------------------|-------------| | Silhouette | 0.084±0.016 | 1.045±0.064 | 0.777±0.116 | 1.808±0.112 | 0.015±0.011 | 1.887±0.118 | 0.032±0.007 | | DBCV | 0.090±0.015 | 1.119±0.073 | 0.897±0.154 | 1.937±0.109 | 0.017±0.012 | 1.934±0.116 | 0.043±0.010 | | No clusters | 0.009±0.004 | 0.453±0.051 | 0.015±0.012 | 1.044±0.104 | 0.005±0.003 | 1.478±0.132 | 0.017±0.011 | Figure 5: Samples generated by ROME and KDE - equivalent to ROME without clustering, decorrelation and normalization - (in pink) contrasted with samples from p (in blue); Aniso. Note that the samples by KDE are more spread out, indicating over-smoothing. ![Image](https://ijcai.org/proceedings/2024/fulltext/525/image_000004_6b9bd9f95f2761e2302a7c6e80d4128c23230d72a4dd4f6f791c4fe9f5c33292.png) However, using clustering is essential, as otherwise there is a risk of over-smoothing, such as in the case of multi-modal distributions with varying densities in each mode (Table 2). Testing variants of ROME on the Aniso distribution (Table 3) demonstrated not only the need for decorrelation, through the use of rotation, but also normalization in the case of distributions with highly correlated features. There, using either of the two clustering methods in combination with normalization and decorrelation (our full proposed method) is better than the two alternatives of omitting only decorrelation or both decorrelation and normalization. In the case of clustering with the silhouette score, the full method is significantly more likely to reproduce the underlying distribution p by a factor of 1.19 (with a statistical significance of 10-50 , see Appendix B) as opposed to omitting only decorrelation, and by 1.30 compared to omitting both decorrelation and normalization. Similar trends can be seen when clustering based on DBCV, with the full method being more likely to reproduce p by a factor of 1.14 and 1.31 respectively. Results on the Aniso distribution further show that KDE on its own is not able to achieve the same results as ROME, but rather it has a tendency to over-smooth (Figure 5). Additionally, the ablation with and without normalization on the Two Moons distribution (Table 4) showed that normalization is necessary to avoid over-smoothing on non-normal distributions. Lastly, investigating the effect of different downstream density estimators, we found that using ROME with fkNN instead of fKDE leads to over-fitting (highest DJS values in Table 5). Meanwhile, ROME with fGMM tends to over-smooth the estimated density in cases where the underlying distribution is not Gaussian (high ̂ W in Table 4). The over-smoothing caused by fGMM is further visualised in Figure 6. In conclusion, our ablation studies confirmed that using fKDE in combination with data clustering, normalization and decorrelation provides the most reliable density estimation for different types of distributions. Figure 6: Samples generated by ROME, and ROME with fGMM as the downstream estimator (in pink) contrasted with samples from p (in blue); Two Moons. Note that the samples from fGMM are more spread out, which clearly displays over-smoothing. ![Image](https://ijcai.org/proceedings/2024/fulltext/525/image_000005_48219c35643d37e61b553f094de52a7a5ca2bc6782f7e568150ae8024332c927.png) ## 6 Conclusion In our comparison against two established and sophisticated density estimators, we observed that ROME achieved consistently good results across all tests, while manifold Parzen windows (MPW) and vine copulas (VC) were susceptible to over-fitting and over-smoothing. For example, while MPW is superior at capturing non-normal distributions, it produces kernels with too small of a bandwidth (hence the over-fitting), which is likely caused by the selected number of nearest neighbors used for the localized kernel estimation being too small. Meanwhile, compared to VC, ROME is numerically more stable and does not hallucinate new modes in the estimated densities (Figure 4). Furthermore, as part of several ablation studies, we found that ROME overcomes the shortcomings of other common density estimators, such as the overfitting exhibited by kNN or the over-smoothing by GMM. In those studies, we additionally demonstrated that our approach of using clustering, decorrelation, and normalization is indispensable for overcoming the deficiencies of KDE. Future work can further improve on our results by investigating the integration of more sophisticated density estimation methods, such as MPW, instead of the kernel density estimator in our proposed approach to enable better performance on non-normal clusters. Overall, by providing a simple way to accurately estimate distributions based on samples, ROME can help in better handling and evaluating probabilistic data as well as enabling more precise probabilistic inference. ## Acknowledgments This research was supported by NWO-NWA project "Acting under uncertainty" (ACT), NWA.1292.19.298. ## Contribution Statement Anna M´esz´aros and Julian F. Schumann equally contributed to the design and implementation of the research, to the analysis of the results and to the writing of the manuscript and should be considered joint first authors. Javier Alonso-Mora provided valuable feedback on the writing of the manuscript. Arkady Zgonnikov and Jens Kober provided valuable feedback at all steps of the project and should be considered joint last authors. ## References 1. [Ankerst et al., 1999] Mihael Ankerst, Markus M Breunig, Hans-Peter Kriegel, and Jörg Sander. Optics: Ordering points to identify the clustering structure. ACM Sigmod Record, 28(2):49-60, 1999. 2. [Bakam and Pommeret, 2023] Yves I Ngounou Bakam and Denys Pommeret. Nonparametric estimation of copulas and copula densities by orthogonal projections. Econo- metrics and Statistics, 2023. 3. [Deisenroth et al., 2020] Marc Peter Deisenroth, A Aldo Faisal, and Cheng Soon Ong. Mathematics for Machine Learning. Cambridge University Press, 2020. 4. [Ester et al., 1996] Martin Ester, Hans-Peter Kriegel, Jörg Sander, and Xiaowei Xu. A density-based algorithm for discovering clusters in large spatial databases with noise. In Proceedings of the Second International Conference on Knowledge Discovery and Data Mining, volume 96, pages 226-231, 1996. 5. [Gao et al., 2022] Jia-Xing Gao, Da-Quan Jiang, and Min- Ping Qian. Adaptive manifold density estimation. Journal of Statistical Computation and Simulation, 92(11):2317- 2331, 2022. 6. [Goodfellow et al., 2020] Ian Goodfellow, Jean Pouget- Abadie, Mehdi Mirza, Bing Xu, David Warde-Farley, Sherjil Ozair, Aaron Courville, and Yoshua Bengio. Gen- erative adversarial networks. Communications of the ACM, 63(11):139-144, 2020. 7. [Heidenreich et al., 2013] Nils-Bastian Heidenreich, Anja Schindler, and Stefan Sperlich. Bandwidth selection for kernel density estimation: a review of fully automatic se- lectors. AStA Advances in Statistical Analysis, 97:403- 433, 2013. 8. [Joe, 2014] Harry Joe. Dependence modeling with copulas. CRC press, 2014. 9. [Kobyzev et al., 2020] Ivan Kobyzev, Simon JD Prince, and Marcus A Brubaker. Normalizing flows: An introduc- tion and review of current methods. IEEE Transactions on Pattern Analysis and Machine Intelligence, 43(11):3964- 3979, 2020. 10. [Kvalseth, 1997] Tarald O Kvalseth. Generalized divergence and gibbs' inequality. In 1997 IEEE International Confer- ence on Systems, Man, and Cybernetics. Computational Cybernetics and Simulation, volume 2, pages 1797-1801. IEEE, 1997. 11. [Liang et al., 2020] Junwei Liang, Lu Jiang, Kevin Mur- phy, Ting Yu, and Alexander Hauptmann. The garden of forking paths: Towards multi-future trajectory predic- tion. In Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition, pages 10508- 10518, 2020. 12. [Lin, 1991] Jianhua Lin. Divergence measures based on the shannon entropy. IEEE Transactions on Information The- ory, 37(1):145-151, 1991. 13. [Loftsgaarden and Quesenberry, 1965] Don O Loftsgaarden and Charles P Quesenberry. A nonparametric estimate of a multivariate density function. The Annals of Mathematical Statistics, 36(3):1049-1051, 1965. 14. [McLachlan and Rathnayake, 2014] Geoffrey J McLachlan and Suren Rathnayake. On the number of components in a gaussian mixture model. Wiley Interdisciplinary Re- views: Data Mining and Knowledge Discovery, 4(5):341- 355, 2014. 15. [Moulavi et al., 2014] Davoud Moulavi, Pablo A Jaskowiak, Ricardo JGB Campello, Arthur Zimek, and Jörg Sander. Density-based clustering validation. In Proceedings of the 2014 SIAM International Conference on Data Mining, pages 839-847. SIAM, 2014. 16. [Mozaffari et al., 2020] Sajjad Mozaffari, Omar Y Al-Jarrah, Mehrdad Dianati, Paul Jennings, and Alexandros Mouza- kitis. Deep learning-based vehicle behavior prediction for autonomous driving applications: A review. IEEE Trans- actions on Intelligent Transportation Systems, 23(1):33- 47, 2020. 17. [Nagler and Czado, 2016] Thomas Nagler and Claudia Czado. Evading the curse of dimensionality in nonpara- metric density estimation with simplified vine copulas. Journal of Multivariate Analysis, 151:69-89, 2016. 18. [Otneim and Tjøstheim, 2017] Håkon Otneim and Dag Tjøstheim. The locally gaussian density estimator for multivariate data. Statistics and Computing, 27:1595- 1616, 2017. 19. [Parzen, 1962] Emanuel Parzen. On estimation of a proba- bility density function and mode. The Annals of Mathe- matical Statistics, 33(3):1065-1076, 1962. 20. [Rasouli, 2020] Amir Rasouli. Deep learning for vision-based prediction: A survey. arXiv preprint arXiv:2007.00095, 2020. 21. [Rousseeuw, 1987] Peter J Rousseeuw. Silhouettes: a graph- ical aid to the interpretation and validation of cluster anal- ysis. Journal of Computational and Applied Mathematics, 20:53-65, 1987. 22. [Silverman, 1998] Bernard W Silverman. Density estimation for statistics and data analysis. Routledge, 1998. 23. [Villani, 2009] Cédric Villani. Optimal transport: Old and new, volume 338. Springer, 2009. 24. [Vincent and Bengio, 2002] Pascal Vincent and Yoshua Ben- gio. Manifold parzen windows. Advances in Neural Infor- mation Processing Systems, 15, 2002. 25. [Wang et al., 2009] Xiaoxia Wang, Peter Tino, Mark A Fardal, Somak Raychaudhury, and Arif Babul. Fast parzen window density estimator. In 2009 International Joint Conference on Neural Networks, pages 3267-3274. IEEE, 2009. 26. [Wei et al., 2020] Ruoqi Wei, Cesar Garcia, Ahmed El- Sayed, Viyaleta Peterson, and Ausif Mahmood. Vari- ations in variational autoencoders-a comparative evalua- tion. IEEE Access, 8:153651-153670, 2020. 27. [Wold et al., 1987] Svante Wold, Kim Esbensen, and Paul Geladi. Principal component analysis. Chemometrics and Intelligent Laboratory Systems, 2(1-3):37-52, 1987. 28. [Xu et al., 2019] Donna Xu, Yaxin Shi, Ivor W Tsang, Yew- Soon Ong, Chen Gong, and Xiaobo Shen. Survey on multi-output learning. IEEE Transactions on Neural Net- works and Learning Systems, 31(7):2409-2429, 2019.