ROME: Robust Multi-Modal Density Estimator

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

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

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.

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]:

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.

• 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 ±.

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.

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.

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.

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.