Bayesian nonparametric methods provide a useful alternative to black box machine learning algorithms, having potential advantages in terms of characterizing uncertainty in inferences and predictions. However, computation can be slow and unwieldy to implement. Hence, it is important to develop simpler and faster Bayesian nonparametric approaches, and hybrid methods that borrow the best of both worlds. For example, if one could use the Bayesian machinery for uncertainty quantification and reduction of mean square errors through shrinkage, while incorporating algorithmic aspects of machine learning approaches, one may be able to engineer a highly effective hybrid. The focus of this article is on proposing such an approach for density estimation, motivated by the successes and limitations of nearest neighbor algorithms and Bayesian mixture models.

Nearest neighbor algorithms are popular due to a combination of simplicity and performance. Given a set of n observations , the density at x is estimated as ˆis the number of neighbors of the distance of x from its kth nearest neighbor in is the volume of the p-dimensional unit ball (Loftsgaarden and Quesenberry, 1965; Mack and Rosenblatt, 1979). Refer to Biau and Devroye (2015) for an overview of related estimators and corresponding theory.

Nearest neighbor density estimators are a type of locally adaptive kernel density estimators. The literature on such methods identifies two broad classes: balloon estimators and sample smoothing estimators; see Scott (2015); Terrell and Scott (1992) for an overview. Balloon estimators characterize the density at a query point x using a bandwidth function h(x); classical examples include the naive k-nearest neighbor density estimator (Loftsgaarden and Quesenberry, 1965) and its modification in Mack and Rosenblatt (1979). More elaborate balloon estimators face challenges in terms of choice of h(x) and obtaining estimators that do not integrate to 1. Sample smoothing estimators make use of n different bandwidths ), one for each sample point , to estimate the density at a query point x globally. By construction, sample smoothing estimators are bona fide density functions integrating to 1. To fit either the balloon or the sample smoothing estimator, one may compute an initial pilot density estimator employing a constant bandwidth and then use this pilot to estimate the bandwidth function (Breiman et al., 1977; Abramson, 1982). Another example of a locally adaptive density estimator is the local likelihood density estimator (Loader, 1996, 2006; Hjort and Jones, 1996), which fits a polynomial model in the neighborhood of a query point x to estimate the density at x, estimating the parameters of the local polynomial by maximizing a penalized local log-likelihood function. The above methods produce a point estimate of the density without uncertainty quantification (UQ).

Alternatively, there is a Bayesian literature on locally adaptive kernel methods, which

express the unknown density as:

which is a mixture of m components, with the hth having probability weight and kernel parameters ; by allowing the location and bandwidth to vary across components, local adaptivity is obtained. A Bayesian specification is completed with prior for the kernel parameters and for the weights. In practice, it is common to rely on an over-fitted mixture model (Rousseau and Mengersen, 2011), which chooses m as a pre-specified finite upper bound on the number of components, and lets

By choosing close to zero, this prior favors effective deletion of redundant components. Also, augmenting with component index , a simple Gibbs sampler can be used for posterior computation, alternating between sampling (i) multinomial conditional posterior, for and (iii) Dirichlet(

Relative to frequentist locally adaptive methods, Bayesian approaches are appealing in automatically providing a characterization of uncertainty in estimation, while having excellent practical performance for a broad variety of density shapes and dimensions. However, implementation typically relies on Markov chain Monte Carlo (MCMC), with the Gibbs sampler sketched above providing an example of a common algorithm used in practice. Unfortunately, current MCMC algorithms for posterior sampling in mixture models tend to face issues with slow mixing, meaning the sampler can take a very large number of iterations to adequately explore different posterior modes and obtain sufficiently accurate posterior summaries.

MCMC inefficiency has motivated a literature on faster approaches, including sequential approximations (Wang and Dunson, 2011; Zhang et al., 2014) and variational Bayes (Blei and Jordan, 2006). These methods are order dependent, tend to converge to local modes, and/or lack theory support. Newton and Zhang (1999); Newton (2002) instead rely on predictive recursion. Such estimators are fast to compute and have theory support, but are also order dependent and do not provide a characterization of uncertainty. Alternatively, one can use a Polya tree as a conjugate prior (Lavine, 1992, 1994), and there is a rich literature on related multiscale and recursive partitioning approaches, such as the optional Polya tree (Wong and Ma, 2010). However, Polya trees have disadvantages in terms of sensitivity to a base partition and a tendency to favor spiky/erratic densities. These disadvantages are

inherited by most of the computationally fast modifications.

This article develops an alternative to current locally adaptive density estimators, obtaining the practical advantages of Bayesian approaches in terms of uncertainty quantification and a tendency to have relatively good performance for a wide variety of true densities, but without the computational disadvantage due to the use of MCMC. This is accomplished with a Nearest Neighbor-Dirichlet Mixture (NN-DM) model. The basic idea is to rely on fast nearest neighbor search algorithms to group the data into local neighborhoods, and then condition on these neighborhoods in defining a Bayesian mixture model-based approach. Section 2 outlines the NN-DM approach and describes implementation details for Gaussian kernels. Section 3 provides some theory support for NN-DM. Section 4 contains simulation experiments comparing NN-DM with a rich variety of competitors in univariate and multivariate examples, including an assessment of UQ performance. Section 5 contains a real data application, and Section 6 a discussion.

2.1 Nearest Neighbor Dirichlet Mixture Framework

Let ) denote a distance metric between data points Euclidean distance is typically chosen. For each denote the jth nearest neighbor to in the data ), such that ), with ties broken by increasing order of indices. The indices on the k nearest neighbors to are denoted as , where by convention we define . Denote the set of data points in the i-th neighborhood by . In implementing the proposed method, we typically let the number of neighbors k vary as a function of n. When necessary we use the notation to express this dependence. However, we routinely drop the n subscript for notational simplicity.

Fixing , we model the density of the data within the i-th neighborhood using

where are parameters specific to neighborhood i that are given a global prior distribution . To combine the )s into a single global f(x), similarly to equations (1)-(2), we let

The key difference relative to standard Bayesian mixture model (1) is that in (4) we include one component for each data sample and assume that only the data in the k-nearest neighborhood of sample i will inform about . In contrast, (1) lacks any sample dependence, and we infer allocation of samples to mixture components in a posterior inference phase.

Given the restriction that only data in the i-th neighborhood inform about pseudo-posterior distribution of given data

This pseudo-posterior is in a simple analytic form if is conjugate to ). The prior can involve unknown parameters and borrows information across neighborhoods; this reduces the large variance problem common to nearest neighbor estimators.

Since the neighborhoods are overlapping, the conditional posterior for not exactly Dirichlet. However, one can define the number of unique points in the i-th neighborhood similar in spirit to the number of points in the h-th cluster in mixture models of the form (1). By convention, we define the point that generated its neighborhood to be a member of that neighborhood. For any other data point to be a unique member of the neighborhood generated by , we require lies in the neighborhood generated by does not lie in the neighborhood of any other . In Section 3.2, we show that the number of unique points defined as above approaches 1 as . This motivates the following Dirichlet pseudo-posterior update for the neighborhood weights

This distribution is equivalent to the conditional posterior on the kernel weights in the Dirichlet mixture of equations (1)-(2), but we use n components and fix the effective number of samples allocated to each component at one. The Dirichlet distribution in (6) is centered on (1/n, . . . , 1/n) with concentration parameter + 1. In practice, is typically chosen to be close to zero, as motivated in Section 1, and we let increase slowly with the sample size.

Based on equations (3)-(6), our nearest neighbor-Dirichlet mixture produces a pseudo-posterior distribution for the unknown density f(x) through simple distributions for the parameters characterizing the density within each neighborhood and for the weights. To generate independent Monte Carlo samples from the pseudo-posterior for f, one can simply draw samples from (5) and (6) independently and plug these samples into the expression for f(x) in (4). Although this is not exactly a coherent fully Bayesian posterior distribution, we claim that it can be used as a practical alternative to such a posterior in practice. This claim is backed up by theoretical arguments, simulation studies and a real data application in the sequel.

2.2 Illustration with Gaussian Kernels

Suppose we have independent and identically distributed (iid) observations density is an unknown density function with respect to the Lebesgue measure on denote the set of all real-valued positive definite matrices. Fix . We will illustrate the method for a general note key changes for the special case p = 1. We proceed by setting ) to be the multivariate Gaussian density . We first compute the neighborhoods corresponding to as in Section 2.1 and place a normal-inverse Wishart prior on ), given by () independently for and Σ; for details about parametrization see Section J of the Appendix.

For p = 1, we have a univariate Gaussian density ) in neighborhood i with normal-inverse gamma priors (2) independently for i = 1IG(= 1, the IW) density simplifies to an IG(density with

Monte Carlo samples from the pseudo-posterior of f(x) given the data obtained using Algorithm 1. The corresponding steps for the univariate case are provided in Section I of the Appendix. Although the pseudo-posterior distribution of f(x) lacks an analytic form, we can obtain a simple form for its pseudo-posterior mean by integrating over the pseudo-posterior distribution of (Recall the definitions of Ψfrom Step 2 of Algorithm 1 and define Λpseudo-posterior mean of f(x) is given by

Here for -dimensional Student’s t-density with degrees of freedom 0, location and scale matrix Λ For the univariate case, we have

Algorithm 1. The pseudo-posterior mean of f(x) is given by

where ) represents the univariate Student’s t-density with degrees of freedom.

2.3 Hyperparameter Choice

The hyperparameters in the prior for the neighborhood-specific parameters need to be chosen carefully – we found results to be sensitive to . If non-informative values are chosen for these key hyperparameters, we tend to inherit typical problems of nearest neighbor estimators including lack of smoothness and high variance. Suppose Σ for denote the i, jth entry of Σ and Ψ, respectively. Then Σ+ 1. Thus borrowing from the univariate case, we set Ψ= 0 for all , which implies that Ψuse leave-one-out cross-validation to select the optimum dimensional data, we recommend fixing which implies a multivariate Cauchy prior predictive density. We choose the leave-one-out log-likelihood as the criterion function for cross-validation, which is closely related to minimizing the Kullback-Leibler divergence between the true and estimated density (Hall, 1987; Bowman, 1984). The explicit expression for the pseudo-posterior mean in equations (7) and (8) makes cross-validation computationally efficient. The description of a fast implementation is provided in Section H of the Appendix.

The proposed method has substantially faster runtime if one bypasses cross-validation and uses a default choice of hyperparameters. In particular, we found the default values to work well across a number of simulation cases, especially when the true density is smooth. Although using cross-validation to estimate can lead to improved performance when the underlying density is spiky, cross-validation provides little to no gains for smooth true densities. Furthermore, with low sample size and increasing number of dimensions, we found this gap to diminish rapidly.

The other key tuning parameter for NN-DM is the number of nearest neighbors The pseudo-posterior mean in (7) reduces to a single . In contrast, = 1 provides a sample smoothing kernel density estimate with a specific bandwidth function (Terrell and Scott, 1992). Therefore, the choice of k can impact the smoothness of the density estimate. To assess the sensitivity of the NN-DM estimate to the choice of k, we investigate how the out-of-sample log-likelihood of a test set changes with respect to k in Section 4.7. These simulations suggest that the proposed method is quite robust to the exact choice of k. In practice with finite samples and small dimensions, we recommend a default choice of = 10 for univariate and multivariate cases, respectively. These values led to good performance across a wide variety of simulation cases as described in Section 4.

3.1 Asymptotic Properties

There is a rich literature on asymptotic properties of the posterior measure for an unknown density under Bayesian models, providing a frequentist justification for Bayesian density estimation; refer, for example to Ghosal et al. (1999), Ghosal et al. (2007). Unfortunately, the tools developed in this literature rely critically on the mathematical properties of fully Bayes posteriors, providing theoretical guarantees for a computationally intractable exact posterior distribution under a Bayesian model. Our focus is instead on providing frequentist asymptotic guarantees for our computationally efficient NN-DM approach, with this task made much more complex by the dependence across neighborhoods induced by the use of a nearest neighbor procedure.

We first focus on proving pointwise consistency of the pseudo-posterior of f(x) induced by (3)-(6) for each , using Gaussian kernels as in Section 2.2. We separately study the mean and variance of the NN-DM pseudo-posterior distribution, first showing that the pseudo-posterior mean in (7) is pointwise consistent and then that the pseudo-posterior variance vanishes asymptotically. The key idea behind our proof is to show that the pseudo-posterior mean is asymptotically close to a kernel density estimator with suitably chosen bandwidth for fixed at a desired rate. The proof then follows from standard arguments leading to consistency of kernel density estimators. The NN-DM pseudo-posterior mean mimics a kernel density estimator only in the asymptotic regime; in finite sample simulation studies (refer to Section 4), NN-DM has much better performance. The detailed proofs of all results in this section are in the Appendix.

Consider independent and identically distributed data from a fixed unknown density with respect to the Lebesgue measure on equipped with the Euclidean metric, inducing the measure to denote the mean of f(x), variance of f(x) and probability of the event ), respectively, under the pseudo-posterior distribution of f(x) implied by equations (3)-(6). We make the following regularity assumptions on

Assumption 3.2 (Bounded gradient)is continuous on [0all and some finite L > 0.

Assumption 3.3 (Bounded sup-norm)

Our asymptotic analysis relies on analyzing the behaviour of the pseudo-posterior updates within each nearest neighborhood. We leverage on key results from Biau and Devroye (2015); Evans et al. (2002) which are based on the assumption that the true density has compact support as in Assumption 3.1. Assumption 3.2 ensures that the kernel density estimator has finite expectation. Versions of this assumption are common in the kernel density literature; for example, refer to Tsybakov (2004). Assumptions 3.1 and 3.3 imply the existence of 0 such that 0 referred to as a positive density condition by Evans (2008); Evans et al. (2002). This is used to establish consistency of the proposed method, justify the choice of the pseudo-posterior distribution of the weights, and obtain our tuning algorithm for the Dirichlet prior concentration parameter . These assumptions are standard in the literature studying frequentist asymptotic properties of nearest neighbor and Bayesian density estimators.

For i = 1, . . . , n, recall the definitions of

where are as in Algorithm 1. Define the bandwidth

matrix

We have suppressed the dependence of for notational convenience. It is immediate that . To prove consistency of the pseudo-posterior mean, we first show that ˆare asymptotically close, that is we show that obtain this result, we approximate using successive applications of the mean value theorem. Finally, we exploit the convergence of ) to the true value to obtain the consistency of ˆ). The proof of convergence of ) is provided in Section G of the Appendix. The precise statement regarding the consistency of the pseudo-posterior mean is given in the following theorem. Let denote the minimum of a and b.

that -probability as

We now look at the pseudo-posterior variance of f(x). We let

For

As 2. Analogous steps to the ones used in the proof of Theorem 3.4 can be used to imply that -probability. Also, as (4using Stirling’s approximation. We now provide an upper bound on the pseudo-posterior variance of f(x) which shows convergence of the pseudo-posterior variance to 0.

be the bandwidth matrix defined in (9). Let . Under Assumptions 3.1-3.3 with as in Theorem 3.4, we have

Refer to Sections B and C in the Appendix for proofs of Theorems 4 and 5, respectively. Pointwise pseudo-posterior consistency follows directly from Theorems 3.4 and 3.5 as shown below.

satisfy Assumptions 3.1-3.3 with as in Theorem 3.4. Fix and define the -ball around denote the probability of the set under the pseudo-posterior distribution of f(x) as induced by equations (3)-(6). Then pr-probability as

0 and consider the . Then by Chebychev’s inequality, we have pr-probability as , using Theorems 3.4 and 3.5.

Our next result focuses on the limiting distribution of g(x), corresponding to a simplified form of the NN-DM which sets the weights to their pseudo-posterior mean,

where we focus on the univariate case with Gaussian kernels. From Section I of the Appendix, the pseudo-posterior distribution of (is given by NIG(are as before and

We show that the limiting distribution of g(x) is a Gaussian distribution with appropriate centering and scaling. This allows interpretation of 100(1 )% pseudo-credible intervals as 100(1)% frequentist confidence intervals on average for large n. Use of simplified versions of the original estimator to derive asymptotic distribution results is standard in the Bayesian nonparametric literature; refer to Chapter 12 of Ghosal and Van der Vaart (2017). Let denote the l-th derivative of

Theorem 3.7. For a fixed denote the simplified NN-DM estimator in (13). Suppose satisfies Assumptions 3.1-3.3 and also satisfies for some finite (9) satisfying

For a proof of Theorem 3.7, we refer the reader to Section D of the Appendix.

3.2 Pseudo-Posterior Distribution of Weights

We investigate the rationale behind the pseudo-posterior update of the weight vector which has a symmetric prior distribution Dirichlet() as motivated in Section 1. As discussed in Section 1, the conditional update for the weights in a finite Bayesian mixture model with m components given the cluster allocation indices is obtained by Dirichlet(is the prior concentration parameter and ) is the number of data points allocated to the h-th cluster. This is not true in our case as the -nearest neighborhoods have considerable overlap between them. Instead, we consider the number of unique data points in each of these neighborhoods.

Define the -nearest neighborhood of to be the set where -th nearest neighbor of in the data , following the notation in Section 2.1. We assume ) is the Euclidean metric from here on, and let denote the distance of -th nearest neighbor in

Let denote the number of unique members in as defined in Section 2.1. Then, we can express

where I(A) is the indicator function of the set

by symmetry. Furthermore, are identically distributed for i = 1, . . . , n. We now state a result which provides a motivation for our choice of the pseudo-posterior update of two sequences of real numbers (), we write some constant

satisfying Assumptions 3.1-3.3. Furthermore, suppose that is as defined in Theorem 3.4.

Proof of Theorem 3.8 is in Section E of the Appendix. The above theorem suggests we asymptotically have only one unique member per neighborhood , namely the point itself generated this neighborhood. This result motivates our choice of the pseudo-posterior update of the weight vector We illustrate uncertainty quantification of the proposed method in finite samples in Section 4.4 with this choice of pseudo-posterior update of the weight vector

3.3 Choice of Dirichlet Prior Parameter

Although Theorem 3.6 implies consistency of the pseudo-posterior for fixed the choice of impacts frequentist coverage of the pseudo-posterior credible intervals as it directly influences the pseudo-posterior variance of f(x) through Theorem 3.5. We now describe a data-dependent method to choose using Bernstein-von Mises results for linear functionals of Bayesian density estimators in the univariate setup (Rivoirard and Rousseau, 2012). The key idea we adopt is to consider a linear functional ) of the density ˜f such as its mean, obtain its Bernstein-von Mises limit distribution using the result in Rivoirard and Rousseau (2012), and equate the variance of this limit distribution with the pseudo-posterior variance of P(f) when f has the NN-DM formulation.

Let ˜satisfy conditions in Rivoirard and Rousseau (2012). Define the linear functional for a generic density ˜independent and identically distributed data from a density supported on [0, 1]. Define be the empirical distribution function of Then under a suitable prior distribution on ˜f, Rivoirard and Rousseau (2012) show that the posterior distribution of is asymptotically Gaussian with mean 0 and variance Ω

As a specific example, consider the case when ˜p(u) = u, corresponding to the linear functional given by . Following Rivoirard and Rousseau (2012), the limiting variance of the posterior distribution of under a suitable prior on ˜f is given by Ωis the population mean and is the population variance of the density . This provides the asymptotic variance of for an appropriate Bayesian density estimator ˜f. Our strategy for finding a value of involves equating the pseudo-posterior variance of when f has the NN-DM formulation, as . This is done to ensure that the same limiting variance as is the proposed estimator and ˜f is a valid Bayesian density estimator according to Rivoirard and Rousseau (2012). When f has the NN-DM formulation with Gaussian kernels, , following from Section 2.2, with the variance of ) provided below.

satisfies Assumptions 3.1-3.3 with p = 1. Let f have the NNDM formulation, and be chosen as in Theorem 3.4. Let Θ =as in Section 2.2. For ¯

The proof of Theorem 3.9 and derivation of (17) are in Section F of the Appendix. Due to the lack of a multivariate analogue of the result discussed in Rivoirard and Rousseau

(2012), a choice of when data are multivariate is not immediate. We observe that the univariate choice of may be written as is large with that end, let be the multivariate bandwidth matrix as defined in Section 3.1 and Σbe the unknown population covariance matrix of . Then, one may potentially extend the findings of (17) to a multivariate extension given by

It is immediate that the choice of in the multivariate case as in (18) reduces to the choice in (17) for is estimated according to Section 2.3 and the underlying population variance is estimated, one can use (17) or its multivariate analogue (18) to select an appropriate value of

4.1 Preliminaries

In this section, we compare the performance of the proposed density estimator with several other standard density estimators through several numerical experiments. We evaluate performance based on the expected distance (Devroye and Gyorfi, 1985). For the pair (is the true data generating density and ˆf is an estimator, the expected distance is defined as . Given a sample size n, we compute an estimate ˆ) in two steps. First, we sample and obtain ˆf based on this sample, and then further sample independent test points and compute

In the second step, to approximate the expectation with respect to , the first step is repeated R times. Letting ˆdenote the estimate for the rth replicate, we compute the final estimate as ˆ. Then, it follows that ˆ, by the law of large numbers. In our experiments, we set We let 0denote the vector with all entries equal to 0 and the vector with all entries equal to 1 in , respectively, for

All simulations were carried out using the R programming language (R Core Team, 2018). For Dirichlet process mixture models, we collect 2, 000 Markov chain Monte Carlo (MCMC) samples after discarding a burn-in of 3, 000 samples using the dirichletprocess package (J. Ross and Markwick, 2019). The default implementation of the Dirichlet process mixture model in p dimensions in the dirichletprocess package uses multivariate Gaussian kernels and has the base measure as NIW) with the Dirichlet concentration parameter having the Gamma(2, 4) prior (West, 1992). For the nearest neighbor-Dirichlet mixture, 1, 000 Monte Carlo samples are taken. For the kernel density estimator, we select the bandwidth by the default plug-in method hpi for univariate cases and Hpi for multivariate cases (Sheather and Jones, 1991; Wand and Jones, 1994) using the package ks (Duong, 2020). We additionally consider the k-nearest neighbor estimator studied in Mack and Rosenblatt (1979), setting the number of neighbors , and the variational Bayes (VB) approximation to Dirichlet process mixture models (Blei and Jordan, 2006). We also compare with the optional Polya tree (OPT) (Wong and Ma, 2010) using the package PTT. For univariate cases, we consider the recursive predictive density estimator (RD) from Hahn et al. (2018), Polya tree mixtures (PTM) using the package DPpackage (Jara et al., 2011), and the sample smoothing kernel density estimator (A-KDE) using the package quantreg. Lastly, we also compare with the local likelihood density estimator (LLDE) using the package locfit for both univariate and multivariate cases. Dirichlet process mixture model hyperparameter values are kept the same in both the MCMC and variational Bayes implementations, with the number of components of the variational family set to 10 for all cases. We denote the nearest neighbor-Dirichlet mixture, Dirichlet process mixture (DPM) implemented with MCMC, kernel density estimator, variational Bayes approximation to the DPM, and k-nearest neigh-

Figure 1: Box plots of ˆ) for the 10 different choices of the true density and different estimators ˆf for univariate data. The box plots for KDE and RD exclude the heavy-tailed cases CA, IE, and SP.

bor density estimator by NN-DM, DP-MC, KDE, DP-VB and KNN, respectively, in tables and figures.

4.2 Univariate Cases

We set denotes the greatest integer less than or equal to . We consider 10 choices of Mildenberger and Weinert, 2012); the specific choices are Cauchy (CA), claw (CW), double exponential (DE), Gaussian (GS), inverse exponential (IE), lognormal (LN), logistic (LO), skewed bimodal (SB), symmetric Pareto (SP), and sawtooth (ST) with default choices of the corresponding parameters. The prior hyperparameter choices for the proposed method are is chosen via the cross-validation method of Section 2.3. Detailed numerical results are deferred to Table 3 in the Appendix. Instead, in Figure 1, we provide a visual summary of the performance of each method under consideration by forming a box plot of the estimated errors of the methods across all the data generating densities. Methods with lower median as indicated by the solid line of the box plot, and smaller overall spread are preferable as they provide higher accuracy and also maintain such accuracy across a collection of true density cases. Results of KNN are omitted in Figure 1 due to much higher values compared to other methods. For the KDE and RD estimator, the plot and the table exclude the results for the heavy-tailed densities CA, IE and SP due to very high

Overall, a major advantage of the proposed method is its versatility among the considered methods. The Bayesian nonparametric methods DP-MC, DP-VB, PTM, OPT, and RD are often close to NN-DM in terms of their performance when the true densities are smooth and do not display locally spiky behavior. However, the NN-DM performs better than other methods in densities where such local behavior is present and performs very close to the best estimator for either the smooth heavy-tailed or thin-tailed densities. The KDE and RD perform well when data are generated from a smooth underlying density. However, there are some cases where the error for KDE and RD is very high. For instance, when n = 500 and is the standard Cauchy (CA) density, the estimated error for the KDE is 38501.85 and the algorithm for the RD estimate did not converge. Both the KDE and RD also perform poorly in very spiky multi-modal densities such as the ST. Compared to the LLDE and the A-KDE, the NN-DM displays similar performance in heavy-tailed and smooth densities when n = 200, with the NN-DM performing better for the spiky densities. However, when n = 500, the NN-DM shows significant improvements over the LLDE and the A-KDE for spiky densities such as the CW and the ST.

In Figure 2, we show the performance of the NN-DM estimator ˆ(with hyperparameters chosen as described earlier) relative to the posterior mean under a DP-MC with default or hand-tuned hyperparameters, when 500 data points are generated from the sawtooth (ST) density. The Dirichlet process mixture with default hyperparameters is unable to detect the multiple spikes, merging adjacent modes to form larger clusters, perhaps due to inadequate mixing of the Markov chain Monte Carlo sampler or to the Gaussian kernels used in the mixture. As a result, we had to hand-tune the hyperparameters for the Dirichlet process mixture to obtain comparable performance with the NN-DM (without hand-tuning). We obtained the best results when changing the hyperparameters of the base measure of the DPMC to NIG(0, 0.01, 1, 1) while keeping the prior on the same as before. This illustrates the deficiency of the DP-MC in estimating densities with spiky local behavior unless we hand-tune the hyperparameters, which requires knowledge of the true density. We also compare the performance of the two methods with a smoother test density in Figure 3, where the data are generated from a skewed bimodal (SB) distribution. Both the estimates are comparable, but the nearest neighbor-Dirichlet mixture provides better uncertainty quantification. Similar results are obtained for n = 1000, and hence are omitted.

4.3 Multivariate Cases

For the multivariate cases, we consider n = 200 and 1000. The number of neighbors is set to k = 10 and the dimension p is chosen from {2, 3, 4, 6}. Recall the definition of from Section 2.2 and let Φ(x) be the cumulative distribution function of the standard Gaussian density. Let . We consider the following cases.

(1)

Figure 2: Plot comparing density estimates for the NN-DM and DP-MC for n = 500 samples generated from the sawtooth (ST) density. Shaded regions correspond to 95% (pseudo) posterior credible intervals. The true density is displayed using dotted lines. The top panel shows the performance of DP-MC with default hyperparameters on the left and with hand-tuned hyperparameters on the right. The bottom panel shows the performance of the NNDM.

Figure 3: Similar to Figure 2, with data of sample size n = 500 generated from the skewed bimodal (SB) density. Left panel shows the NN-DM fit and the right panel shows the DP-MC fit.

(2) is the diagonal matrix with diagonal entries . We choose and the skewness parameter vector (3) ) is the density of the p-dimensional multivariate Student’s t-distribution as in Section 2.2. We set (4) Mixture of multivariate skew t-distributions (MST): We consider a two component mixture of multivariate skew t-distribution (Azzalini, 2005) given by 0) is the skew t-density with parameters defined as before and the same as in the first case. (5) (6) where denote the density and distribution function of the univariate gamma distribution with shape parameter and rate parameter , respectively, for j = 1, . . . , p and Γ) is as described in Song (2000). This is a Gaussian copula based construction of the multivariate gamma distribution. We set

The hyperparameters for the nearest neighbor-Dirichlet mixture are chosen as 0, where the optimal is chosen via cross-validation as described in Section 2.3. Default hyperparameters as described in Section 4.1 are chosen for the MCMC and VB implementations of the DPM.

Similar to the univariate case, we defer the numerical results to Table 4 in the Appendix and in Figure 4 display a visual summary consisting of box plot of errors over the densities considered. The proposed method is very robust against a wide selection of true distributions, with its error scaling nicely with the dimension. The KDE shows a noticeably sharp decline in performance - when the dimension is changed from 2 to 6, the average increase in error is by factors of about 5 and 7 for sample sizes 200 and 1000, respectively. This is possibly due to lack of adaptive density estimation in higher dimensions using a single bandwidth matrix, since data in become increasingly sparse with increasing p. As in the univariate case, we had to exclude the MVC density for the KDE due to the algorithm not converging. The performances of NN-DM, DP-MC, and DP-VB are quite competitive across densities, with NN-DM faring better than the DP-VB when estimating densities such as the MVC and the MVG. Furthermore, the NN-DM is hit the least significantly by the curse of dimensionality out of the three. This is particularly prominent when n = 200 and p = 6 for the DP-MC when the true density is either MG or MST, and for the DP-VB when the true density is MVC. It is also important to keep in mind that the NN-DM provides similar results compared to the DP-MC while being at least an order of magnitude faster, as illustrated in Section 4.6. The performance of the OPT is hit quite significantly as the number of dimensions increases, along with the algorithm not converging for p = 6. The

Figure 4: Box plots of ˆ) for the 6 different choices of the true density and different estimators ˆf for multivariate data. The box plots for KDE and LLDE exclude the MVC density. The box plots for p = 6 exclude results from OPT.

LLDE provides competitive results with the NN-DM in lower dimensions. However, in higher dimensions, the LLDE often does not converge, indicating lack of stability of the algorithm. We reported the average of the replicates for which the algorithm did converge. The results suggest that the performance of the LLDE is also affected quite drastically with increasing dimensions. When compared across all data generating cases considering the variation in densities, dimensions and sample sizes, the proposed method is seen to be more versatile than its competitors.

4.4 Accuracy of Uncertainty Quantification

In this section, we assess frequentist coverage of 95% pseudo-posterior credible intervals for the NN-DM and compare with coverage based on the 95% posterior credible intervals obtained from DP-MC and DP-VB. Ghosal and Van der Vaart (2017) recommend investigating the frequentist coverage of Bayesian credible intervals. We do not include frequentist coverage for Polya tree mixtures (PTMs) and the optional Polya tree (OPT) due to the lack of available code. We consider the cases in our experiments with sample size n = 500. For each choice of density = 200 test points generated from the density . With these fixed test points, we generate n = 500 data points in our sample for = 200 times and check the coverage of posterior/pseudo-posterior credible intervals obtained from the three methods. We implement the DP-MC with base measure NIW) and a Gamma(2, 4) prior on the concentration parameter as in West (1992). These choices of hyperparameters were seen to give better frequentist coverage results than using the default values used in Sections 4.2 and 4.3. Same choices of hyperparameters are maintained for DP-VB. For the NN-DM we take k = 8 in the univariate case, k = 5 in the bivariate case, as in Section 3.3 and other hyperparameters chosen as before. In Table 1 and in Table 2, we report the average coverage probability and average length of the (pseudo) credible intervals across all the points in the test data for the univariate and bivariate cases, respectively.

For univariate densities, both the DP-MC and DP-VB display severe under-coverage. In most of the cases, the DP-VB and NN-DM have similar width of (pseudo) credible intervals but the DP-VB displays dramatically lower coverage than the NN-DM. The under-coverage displayed by the DP-MC may be due to MCMC mixing issues. The NN-DM shows near nominal coverage in the smooth Gaussian (GS) and lognormal (LN) densities, while also attaining near nominal coverage in the skewed bimodal (SB), claw (CW) and sawtooth (ST) densities which are multi-modal. The shortcomings of DP-MC and DP-VB are especially noticeable when dealing with spiky densities such as the claw or sawtooth. For bivariate cases considered in Table 2 we see a similar trend; the NN-DM method provides uniformly better uncertainty quantification across all the densities considered. It is clear that in terms

Table 1: Comparison of the frequentist coverage of 95% (pseudo) posterior credible intervals of the nearest neighbor-Dirichlet mixture and the MCMC and variational implementations of the Dirichlet process mixture for univariate data. Average length of the intervals are also provided for each case within parentheses. Number of replications and sample size are = 500, respectively.

Table 2: Comparison of the frequentist coverage of 95% (pseudo) posterior credible intervals of the nearest neighbor-Dirichlet mixture and the MCMC and variational implementations of the Dirichlet process mixture for bivariate data. Average length of the intervals are also provided for each case within parentheses. Number of replications and sample size are = 500, respectively.

of frequentist uncertainty quantification, the NN-DM displays vastly superior coverage to the DP-MC and the DP-VB without inflating the interval width.

4.5 Comparison for high dimensional data

In addition to the above experiments, we performed a simulation experiment for high-dimensional data. Specifically, we set p = 50 and consider the same set of true densities in Section 4.3. We compared results from the proposed NN-DM method and the DP-VB. Due to severe computational time, we did not consider the DP-MC in this scenario. We

Figure 5: Out-of-sample log-likelihood of NN-DM and DP-VB on a test set of 500 points for 6 different multivariate densities considered in Section 4.3, for n = 1000 and p = 50.

also tried optional Polya trees (Wong and Ma, 2010) using the PTT package; however, the current implementation of the method breaks down in this high-dimensional setup. Due to numerical instability in computing the error in higher dimensions, we evaluate the methods in terms of their out-of-sample log-likelihood (OOSLL) instead (Gneiting and Raftery, 2007), on a test set of 500 data points. We report the average OOSLL over 30 replications in Figure 5. The results indicate that both methods perform very similarly in terms of out-of-sample fit to the data, with the NN-DM outperforming the DP-VB when the true density is MVC. We also observed that for this experiment, the NN-DM methods with default choice of hyperparameters and with cross-validated choice of Ψhave almost identical performance. For the NN-DM, we set k = 12 after carrying out a sensitivity analysis on k by considering k = 5, 7, 10, 15, and 20. The best results for the NN-DM were obtained for with negligible difference in out-of-sample log-likelihoods between these three choices, with k = 12 performing the best.

4.6 Runtime Comparison

The results in the previous subsection suggest that the DP-VB has dramatic under-coverage with respect to the nominal frequentist coverage. Hence, the VB approach is not useful for uncertainty quantification. For this reason, we focus our comparison of runtimes on the NN-DM and the DP-MC, noting that VB algorithms are faster than either of these approaches.

With n data points in p dimensions, the initial nearest neighbor allocation into n neighborhoods can be carried out in O(n log n) steps (Vaidya, 1986; Ma and Li, 2019). Once the neighborhoods are determined with points in each neighborhood, obtaining the neighborhood specific empirical means and covariance matrices has complexity. Obtaining the pseudo-posterior mean (7) then requires inversion of matrices to evaluate the multivariate t-density, with a runtime of ). Therefore, the total runtime to obtain the pseudo-posterior mean is of the order are interested in uncertainty quantification, we require Monte Carlo samples of the NN-DM, which are independently drawn from its pseudo-posterior. This involves sampling the Dirichlet weights, the neighborhood specific unknown mean and covariance matrix parameters of the Gaussian kernel, and evaluating a Gaussian density for each neighborhood, as outlined in Algorithm 1. To obtain M Monte Carlo samples, the combined complexity of this step is thus ). Overall the runtime complexity to obtain NN-DM samples is therefore ). For high dimensional scenarios, this runtime can be greatly improved by using a low rank matrix factorization of both the neighborhood specific empirical covariance matrix and the sampled covariance matrix parameter to make matrix inversion more efficient (Golub and van Loan, 1996). We now provide a detailed simulation study of runtimes of the proposed method.

In our experiments, we focus on p = 1 and p = 4. The runtime for NN-DM consists of the time to estimate by cross-validation as in Section 2.3 and then drawing samples from its pseudo-posterior. For both dimensions, the sample size is varied from n = 200 to n = 1500 in increments of 100. Data are generated from the standard Gaussian density (GS) for p = 1 and from a mixture of skew t-distributions with the parameters as described for the case MST in Section 4.3 for p = 4. For p = 1, we evaluate the two methods at 500 test points, while for p = 4 we evaluate the methods at 200 test points. The hyperparameters are kept the same as in Sections 4.2 and 4.3. We took 1000 Monte Carlo samples for the NN-DM and 2500 MCMC samples for the DP-MC with a burn-in of 1500 samples. The simulations were carried out on an i7-8700K processor with 16 gigabytes of memory.

In the top panel of Figure 6, we plot the average of the logarithm of the run times of each approach for 10 independent replications. The corresponding average error of the two methods is also included in the bottom panel of Figure 6. The NN-DM is at least an order of magnitude faster than DP-MC. The time saved becomes more pronounced in the multivariate case, where for sample size 1500 the NN-DM is 15 times faster. The gain in computing time does not come at the cost of accuracy as can be seen from the right panel; the proposed method maintains the same order of error as the DP-MC in the univariate case and often outperforms the DP-MC in the multivariate case. We did not implement the Monte Carlo sampler for the proposed algorithm in parallel, but such a modification would substantially improve runtime. Bypassing cross-validation and choosing default hyperparameters instead as outlined in Section 2.3, NN-DM took 7.7 seconds and 28.4 seconds when p = 1 and p = 4, respectively, with sample size n = 1500. In the same

Figure 6: Runtime comparison of DP-MC and NN-DM in univariate case and for 4-dimensional data. Top panel shows runtimes in log scale whereas bottom panel shows corresponding error. Sample size n is varied from 200 to 1500 in increments of 100.

scenario, DP-MC took 291.3 seconds and 1504.4 seconds for p = 1 and p = 4, respectively. Thus the NN-DM with default hyperparameters is about 38 times faster when p = 1 and about 53 times faster when p = 4.

4.7 Sensitivity to the choice of k

In this subsection, we investigate the role of in finite samples for the proposed method. We consider n = 200 samples from the SP density in the univariate case and the MG density in the bivariate case. In each case, we fix a test set of = 500 points, and evaluate the out-of-sample log-likelihood (OOSLL) of the test points for 20 different integer values of k ranging from 2 to 50. Finally, we report results averaged from 10 independent replicates of this setup. We note that for each considered value of k, the parameter estimated using leave-one-out cross-validation. Figure 7 shows how the OOSLL averaged over replicates changes as a function of k for each density considered. The original OOSLL values of the test data points were scaled by the number of test points = 500 for better representability.

For the univariate SP density, the optimal value of k which maximizes the average OOSLL is = 9. This is close to the choice of k = 6 as taken in Section 4.2. For the bivariate MG density, we observe that the choice of k maximizing the OOSLL is = 12, which is also close to the choice of k = 10 as taken in Section 4.3. For both the univariate and the bivariate

Figure 7: Average out-of-sample log-likelihood of 500 test points for the NN-DM as a function of k for one-dimensional and two-dimensional data. Number of samples and number of replications are n = 200 and R = 10, respectively.

case, the out-of-sample log-likelihood of the test set shows little variation with changing k. This indicates that the estimates obtained from the proposed method are quite robust to the particular choice of k.

We apply the proposed density estimator to binary classification. Consider data D = -dimensional feature vectors and are binary class labels. To predict the probability that = 1 for a test point

Bayes rule:

where ˜) is the feature density at ) is the marginal probability of class j, for j = 0, 1. Based on test data, we let = 1), and use either the NN-DM pseudo-posterior mean ˆDP-MC posterior mean ˆ), or the DP-VB posterior mean ˆ) for estimating the within class densities. We omit the kernel density estimator as to the best of our knowledge, no routine R implementation is available for data having more than 6 dimensions. We compare the resulting classification performances in terms of sensitivity, specificity and probabilistic calibration.

The high time resolution universe survey data (Keith et al., 2010) contain information on sampled pulsar stars. Pulsar stars are a type of neutron stars and their radio emissions are detectable from the Earth. These stars have gained considerable interest from the scientific community due to their several applications (Lorimer and Kramer, 2012). The data are

Figure 8: Sensitivity and specificity of the NN-DM, DP-MC, and DP-VB for the universe survey data.

publicly available from the University of California at Irvine machine learning repository. Stars are classified into pulsar and non-pulsar groups according to 8 attributes (Lyon, 2016). There are a total of 17898 instances of stars, among which 1639 are classified as pulsar stars. We create a test data set of 200 stars, among which 23 are pulsar stars. The training size is then varied from 300 to 1800 in increments of 300, each time adding 300 training points by randomly sampling from the entire data leaving out the initial test set. In Figure 8, we plot the sensitivity and specificity of the three methods in consideration. All the methods exhibit similar sensitivity across various training sizes; the DP-MC has marginally better specificity for training sizes 1200 and 1500, while the NN-DM has better specificity for training sizes 300 and 600. Both the NN-DM and the DP-MC exhibit higher specificity and sensitivity than the DP-VB across all training sample sizes considered. We also compare the methods using the Brier score, a proper scoring rule (Gneiting and Raftery,

2007) for probabilistic classification. Suppose for test points and the ith Monte Carlo sample, denotes the sampled 1 probability vector for a generic method. We compute the normalized Brier score for the ith sample as (1is the vector of class labels in the test set. Then with T samples of , we compute the mean Brier score for the three methods considered. The mean Brier score for each training size is shown in the right panel of Figure 9, which naturally shows a declining trend with increasing training size. There is little to choose between the three classifiers in terms of mean Brier score; the proposed method fairs equally well in terms of calibration of estimated test set probabilities with the MCMC implementation of the Dirichlet process. In the left panel of Figure 9, the receiver operating characteristic curve of the methods is shown for 1800 training samples. The area under the curve (AUC) for the NN-DM, the DP-MC and the DP-VB are 0.96, 0.95 and 0.96, respectively. For 1800 training samples, the computation time for the proposed method is about 13 minutes while for the DP-MC it is approximately 5 hours.

Figure 9: Left plot shows the receiver operating characteristic curve of the NN-DM, DP-MC, and DP-VB with 1800 training samples. Area under the curve is abbreviated as AUC. Right plot shows normalized Brier scores for the methods with varying training sample size.

Figure 10: Left and right plot show the out-of-sample log-likelihoods of NN-DM, DP-MC, and DP-VB for the two different star types.

Hence, the proposed method is much faster, even without exploiting parallel computation. We also fitted the proposed method using the training set of all 17698 points; DP-MC was too slow in this case. The sensitivity and specificity of the proposed method increased to 0.99 and 0.91, respectively. We additionally evaluated the methods in terms of the out-of-sample log-likelihood. The results are displayed in Figure 10. While the methods perform comparably in terms of their classification performance, NN-DM achieves a better fit overall, especially for the significantly less prevalent pulsar star type.

The proposed nearest neighbor-Dirichlet mixture provides a useful alternative to Bayesian density estimation based on Dirichlet mixtures with much faster computational speed and stability in avoiding MCMC. MCMC can have very poor performance in mixture models and other multimodal cases, due to difficulty in mixing, and hence can lead to posterior inferences that are unreliable. There is a recent literature attempting to scale up MCMCbased analyses in model-based clustering contexts including for Dirichlet process mixtures; refer, for example to Song et al. (2020) and Ni et al. (2020). However, these approaches are complex to implement and are primarily focused on the problem of clustering, while we are instead focused on flexible modeling of unknown densities.

The main conceptual disadvantage of the proposed approach is the lack of a coherent Bayesian posterior updating rule. However, we have shown that nonetheless the resulting pseudo-posterior can have appealing behavior in terms of frequentist asymptotic properties, finite sample performance, and accuracy in uncertainty quantification. In addition, it is important to keep in mind that Bayesian kernel mixtures have key disadvantages that are difficult to remove within a fully coherent Bayesian modeling framework. These include a strong sensitivity to the choice of kernel and prior on the weights on these kernels; refer, for example to Miller and Dunson (2019).

There are several important next steps. The first is to develop fast and robust algorithms for using the nearest neighbor-Dirichlet mixture not just for density estimation but also as a component of more complex hierarchical models. For example, one may want to model the residual density in regression nonparametrically or treat a random effects distribution as unknown. In such settings, one can potentially update other parameters within a Bayesian model using Markov chain Monte Carlo, while using algorithms related to those proposed in this article to update the nonparametric part conditionally on these other parameters.

R package NNDM available at https://github.com/shounakchattopadhyay/NN-DM was used for the numerical experiments. This research was partially supported by grants R01ES027498 and R01ES028804 of the United States National Institutes of Health and grant N00014-16-1-2147 of the Office of Naval Research.

We first introduce some notation with accompanying technical details which will be used hereafter. We define the Frobenius norm of the with real-valued entries by and denote its determinant by |A|. We observe that, for a vector is the Euclidean norm of a. For two symmetric , we say that is positive semi-definite, that is . For a symmetric matrix , let the eigenvalues of be denoted by ), arranged such that , then it follows by the min-max theorem (Teschl, 2009) that, for each ). In particular, we have

Now consider independent and identically distributed data supported on the interval [0and satisfying Assumptions 3.1-3.3 as mentioned in Section 3.1. Let () and suppose induces the measure on the Borel , denoted by ). We form the k-nearest neighborhood of using the Euclidean norm for i = 1, . . . , n. For a generic -th nearest neighbor in and let be the distance between . Define the ball and the probability of the ball Let denote the sample points which fall in . Then, we define the neighborhood specific empirical mean and covariance matrix as ¯and , respectively. Note that the random vector () is identically distributed for are independent and identically distributed. Thus we only consider the case i = 1 from here on. For sake of brevity, denote by let k depend on n and express this dependence as when required. However, we routinely drop this dependence for notational simplicity.

Conditional on 0, following Mack and Rosenblatt (1979) the conditional joint density of

where ) denotes the indicator function of the event ). Thus conditional on , the random variables are independent and identically distributed, and independent of Q. Also, Mack and Rosenblatt (1979) states that under Assumptions 3.1-3.3 which imply is bounded and continuous on [0have

as is the volume of the unit ball in . Let the function is a non-negative integer. This function can be identified with ) in equation (11) of Mack and Rosenblatt (1979). In the following propositions we will require the expected values of ) for different choices of end, we shall repeatedly make use of the equation (12) from Mack and Rosenblatt (1979) adapted to our setting:

Suppose are independent and identically distributed random variables generated from the density supported on [0satisfying Assumptions 3.1-3.3. For i = 1, . . . , n, recall the definitions of from equation (7):

We want to show that ˆ-probability as ) is as described in (7). We first prove two propositions involving successive mean value theorem type approximations to ˆ), which will imply the final result. We now state the two propositions, with accompanying proofs, before stating the final theorem.

Proof. Since the (Λare identically distributed and (are identically distributed, we have tivariate mean value theorem now implies that

If we let 1)following the choice of Ψfrom Section 2.3, then it is clear that Λ. Therefore, we have Straightforward calculations show that where . Using Theorem 2.4 from Biau and Devroye (2015) for (22) for p = 1, one gets

for an appropriate constant 0. Thus, we have

for sufficiently large n. This implies that

We also have . Finally, simple calculations yield that

where 0 satisfies . Plugging all these back in (23), we obtain a finite constant 0 such that

which goes to 0 as , completing the proof.

We now provide the second mean value theorem type approximation which approximates the random bandwidth matrix Λ

. Also, let . Then, we have

Proof. Using the identically distributed properties of (Λ, we obtain ). Using the multivariate mean value

theorem, we obtain that

where for some Σin the convex hull of Λ, we immediately have Σas well. Using the definitions of Λwe have

Since ¯, we get for sufficiently large n the following:

using (24) and . Taking partial derivatives of logwith respect to Σ evaluated at Σand taking Frobenius norm of both sides, we obtain

for sufficiently large n. We now observe that

where 0. Note that

as 1. This immediately implies that for sufficiently large n. Plugging all these back in equation (27), we obtain for sufficiently large 0 such that

which goes to 0 as , proving the proposition.

We now prove Theorem 3.4.

triangle inequality. Using Propositions B.1 and B.2, we obtain that 0 as . From Section G of the Appendix, we obtain -probability. This immediately implies that given the conditions on and for any ˆ-probability.

) and suppose Then, we have Dirichlet(1+ 1) given . We begin with the identity

For i = 1, . . . , n, we have

where

where the first equality is obtained using last equality is obtained using equation (32). For

We now analyze the second term on the right hand side of (31). Recall that conditionally independent of given the data following from the nearest neighborDirichlet mixture framework as discussed in Section 2.1. Let Σdenote the pseudo-posterior covariance matrix of Dirichlet(+1), standard results yield that Σidentity matrix, and . Then, we have

Using the expression for Σalong with equation (34), we obtain,

where ¯. We now have

where the last inequality is obtained using equation (32). Using as before, we have

Combining equations (33) and (36) and putting the results back in equation (31), we have the inequality. If we let we immediately obtain that var-probability.

Proof. We have iid data satisfying Assumptions 3.1-3.3 for p = 1. The simplified NN-DM density estimator is given by

where (2). The pseudo-posterior distribution of g(x) given is induced through the pseudo-posterior distributions of (. The pseudo- posterior mean is of the form ˆ1)

for from Section B of the Appendix, where

We want to investigate the asymptotic distribution of . For that, let us start with the asymptotic distribution of ), which can be expressed as . Using Lyapunov’s central limit theorem and denoting convergence in distribution under

for some By standard calculations, we have

For r = 3,

It is straightforward to see that (1) for any Lyapunov’s condition is satisfied as the ratio in this case satisfies . Additionally, 0. So by a combination of Lyapunov’s central limit theorem and Slutsky’s theorem, we have

since. From the calculations in Section G of the Appendix, we can expand the Taylor series to two more terms to obtain

since ). The pseudo-posterior variance of g(x) is given by

with being the normalizing constant of the Student’s t-density with degrees of freedom d > 0. Using Stirling’s approximation, ∆This immediately implies

where

Using the techniques of Section B of the Appendix, it can be shown that as . Therefore, we have

as . A simple application of Chebychev’s inequality implies (in -probability as Combining this with (37) and (40) and using Slutsky’s theorem, we obtain the desired result. This result means that pseudo-credible intervals can be considered frequentist confidence intervals on average asymptotically.

E.1 A property of the -nearest neighbor distance

Suppose a density on satisfying Assumptions 3.1-3.3. Denote the induced probability measure for sake of convenience. We define the smoothed k-nearest neighborhood of the Euclidean distance between -nearest neighbor of amongst the data points . It is immediate that are identically distributed from symmetry. Suppose and define the quasi-neighborhood ˜, where the random variables have been replaced by

The positive density condition on obtained from Assumptions 3.1 and 3.3 (Evans et al., 2002; Evans, 2008) ensures the existence of 0 such that for all 0

for all

We first state a Lemma proving some important properties of . Recall that two sequences () are said to be asymptotically equivalent if denoted by

Proof. (i) Note that for sufficiently large n. From Lemma 4.1 of Evans et al. (2002) we have

(ii) For for a sequence Θ0. This ensures that

(iii) Since , a direct application of the first Borel-Cantelli lemma proves the statement.

We now use the above Lemma to prove Theorem 3.8. The key idea is to leverage the fact that with probability 1 for all but finite n.

E.2 Number of unique points in each neighborhood

We now prove Theorem 3.8.

Proof. Using (iii) from Lemma E.1, for i = 1, . . . , n, we have an integer such that for all ) = 1. However, since are identically distributed,

, say. Thus, for all ) = 1 for all immediately implies that

] = 1 for all , which shows ] = 1. Therefore, we have

where for all . Given the conditions on k, it follows that as

for all 0. Therefore, we have

as . This proves the result.

F.1 Proof of Theorem 3.9

. Then, we have Θ = (Dirichlet(start out by observing that

Let n be sufficiently large so that

Since the pseudo-posterior covariance matrix of obtain

where ¯. Putting these back in equation (44) we get that

Combining the results of equations (43) and (45), putting them back in equation (42), and multiplying both sides by n, we get the result.

F.2 Choice of

Suppose is the variance of the underlying true density satisfying Assumptions 3.1-3.3. Let . We start out by observing that

by the triangle inequality and the fact that identically distributed and (are identically distributed, we have

But ) by the triangle inequality, from which it follows that ) are identically distributed for i = 1, . . . , n. Thus, we have

using (25). Since -probability as by the weak law of large numbers, we get that -probability as as well. Equating the pseudo-posterior variance of Θ from Theorem 3.9 with , we get after some rearranging,

As , since each -probability for i = 1, . . . , n using equation (29) with is well approximated by () in the sense that (-probability as . In particular, we have ¯probability. Combining these with the fact that -probability as obtain,

Since 1/n can be asymptotically neglected in comparison to ¯owing to the fact that ) implies the choice of as described in equation (17).

Define the standard multivariate t-density with d > 0 degrees of freedom to be as defined in Section 3.1 is diagonal, it immediately follows that . The following lemma proves the consistency of any such generic kernel density estimator with t kernel depending on n, say

where the bandwidth , with independent and identically distributed data satisfying Assumptions 3.1-3.3.

is a sequence satisfying Let -probability

Proof. It is enough to show that Let us start first with

using the mean value theorem and Polya’s theorem (P´olya, 1920) along with Assumption 3.2 to bound , this implies that

The variance may be dealt with in a similar manner. Following the same steps as before we get

which shows that the variance goes to 0 as

For the nearest neighbor-Dirichlet mixture, recall from Section 3.1 of the main document, where 1)(. Here, the bandwidth G.1 then shows that ) converges to -probability as

H.1 Algorithm for leave-one-out cross-validation

Consider independent and identically distributed data the nearest neighbor-Dirichlet mixture formulation. In the univariate setting, the prior for each of the neighborhood specific parameters is The equivalent prior in the general multivariate setting following Sections 2.2 and 2.3 is (pseudo-posterior mean in equations (7) and (8) to compute leave-one-out log-likelihoods ) for different choices of the hyperparameter maximize this criteria. The details of the computation of ) for a fixed are provided in Algorithm 2.

H.2 Fast Implementation of cross-validation

In Algorithm 2, the nearest neighborhood specification for each is different for i = 1, . . . , n. However, we bypass this computation by initially forming a neighborhood of size (k + 1) for each data point using the entire data and storing the respective neighborhood means and covariance matrices. Suppose for , the indices of the (k + 1)-nearest neighbors are given by ranged in increasing order according to their distance from the neighborhood mean and the neighborhood covariance ma- trix . Then, to form a k-nearest neighborhood for the new data , a single pass over the initial neighborhoods is sufficient to update the new neighborhood means and covariance matrices. Below, we describe the update for the neighborhood means and covariance matrices considering the data

Suppose we have independent and identically distributed observations from the density . In the nearest neighbor-Dirichlet mixture framework for univariate data with the Gaussian kernel ), neighborhood specific parameters independently for i = 1, . . . , n. Monte Carlo samples for the estimated density at any point x can be generated following the steps of Algorithm 3.

The parametrization of the inverse Wishart density defined on the set of all with real entries used in this article is given as follows. Suppose 1 and Ψ is a positive definite matrix. If Σ Ψ), then Σ has the following density function:

where Γ) is the multivariate gamma function defined by

for 2 and the function etr (A) = exp {tr(A)} for a square matrix A. When p = 1, the IWΨ) density is the same as the IG(2) density, where

Table 3: Comparison of the methods in terms of error in the univariate case. Number of test points and replications considered are = 20, respectively.

Abramson, I. S. (1982). On bandwidth variation in kernel estimates-a square root law. The Annals of Statistics, pages 1217–1223.

Azzalini, A. (2005). The skew-normal distribution and related multivariate families. Scandinavian Journal of Statistics, 32(2):159–188.

Biau, G. and Devroye, L. (2015). Lectures on the Nearest Neighbor Method. Springer.

Blei, D. M. and Jordan, M. I. (2006). Variational inference for Dirichlet process mixtures. Bayesian Analysis, 1(1):121–143.

Bowman, A. W. (1984). An alternative method of cross-validation for the smoothing of density estimates. Biometrika, 71(2):353–360.

Breiman, L., Meisel, W., and Purcell, E. (1977). Variable kernel estimates of multivariate densities. Technometrics, 19(2):135–144.

Devroye, L. and Gyorfi, L. (1985). Series in Probability and Statistics.

Duong, T. (2020). ks: Kernel Smoothing. R package version 1.11.7.

Evans, D. (2008). A law of large numbers for nearest neighbour statistics. Proceedings of the Royal Society A: Mathematical, Physical and Engineering Sciences, 464(2100):3175–3192.

Evans, D., Jones, A. J., and Schmidt, W. M. (2002). Asymptotic moments of near–neighbour distance distributions. Proceedings of the Royal Society of London. Series A: Mathematical, Physical and Engineering Sciences, 458(2028):2839–2849.

Ghosal, S., Ghosh, J. K., Ramamoorthi, R., et al. (1999). Posterior consistency of Dirichlet mixtures in density estimation. The Annals of Statistics, 27(1):143–158.

Ghosal, S. and Van der Vaart, A. (2017). Fundamentals of nonparametric Bayesian inference, volume 44. Cambridge University Press.

Ghosal, S., Van Der Vaart, A., et al. (2007). Posterior convergence rates of Dirichlet mixtures at smooth densities. The Annals of Statistics, 35(2):697–723.

Gneiting, T. and Raftery, A. E. (2007). Strictly proper scoring rules, prediction, and esti- mation. Journal of the American Statistical Association, 102(477):359–378.

Golub, G. H. and van Loan, C. F. (1996). Matrix Computations. John Hopkins University Press, 3 edition.

Hahn, P. R., Martin, R., and Walker, S. G. (2018). On recursive bayesian predictive distri- butions. Journal of the American Statistical Association, 113(523):1085–1093.

Hall, P. (1987). On Kullback-Leibler loss and density estimation. The Annals of Statistics, 15(4):1491–1519.

Hjort, N. L. and Jones, M. C. (1996). Locally parametric nonparametric density estimation. The Annals of Statistics, pages 1619–1647.

J. Ross, G. and Markwick, D. (2019). dirichletprocess: Build Dirichlet Process Objects for Bayesian Modelling. R package version 0.3.1.

Jara, A., Hanson, T., Quintana, F., M¨uller, P., and Rosner, G. (2011). DPpackage: Bayesian semi- and nonparametric modeling in R. Journal of Statistical Software, 40(5):1–30.

Keith, M., Jameson, A., Van Straten, W., Bailes, M., Johnston, S., Kramer, M., Possenti, A., Bates, S., Bhat, N., Burgay, M., et al. (2010). The High Time Resolution Universe Pulsar Survey I, System configuration and initial discoveries. Monthly Notices of the Royal Astronomical Society, 409(2):619–627.

Lavine, M. (1992). Some aspects of Polya tree distributions for statistical modelling. Annals of Statistics, 20(3):1222–1235.

Lavine, M. (1994). More aspects of Polya tree distributions for statistical modelling. The Annals of Statistics, 22(3):1161–1176.

Loader, C. (2006). Local regression and likelihood. Springer Science & Business Media.

Loader, C. R. (1996). Local likelihood density estimation. The Annals of Statistics, 24(4):1602–1618.

Loftsgaarden, D. O. and Quesenberry, C. P. (1965). A nonparametric estimate of a multi- variate density function. The Annals of Mathematical Statistics, 36(3):1049–1051.

Lorimer, D. R. and Kramer, M. (2012). Handbook of Pulsar Astronomy.

Lyon, R. J. (2016). Why are pulsars hard to find? PhD thesis, The University of Manchester (United Kingdom).

Ma, H. and Li, J. (2019). A true O(n log n) algorithm for the all-k-nearest-neighbors problem. In International Conference on Combinatorial Optimization and Applications, pages 362– 374. Springer.

Mack, Y. and Rosenblatt, M. (1979). Multivariate k-nearest neighbor density estimates. Journal of Multivariate Analysis, 9(1):1–15.

Mildenberger, T. and Weinert, H. (2012). The benchden package: Benchmark densities for nonparametric density estimation. Journal of Statistical Software, 46(14):1–14.

Miller, J. W. and Dunson, D. B. (2019). Robust Bayesian inference via coarsening. Journal of the American Statistical Association, 114(527):1113–1125.

Newton, M. A. (2002). On a nonparametric recursive estimator of the mixing distribution. , 64(2):306–322.

Newton, M. A. and Zhang, Y. (1999). A recursive algorithm for nonparametric analysis with missing data. Biometrika, 86(1):15–26.

Ni, Y., Ji, Y., and M¨uller, P. (2020). Consensus monte carlo for random subsets using shared anchors. Journal of Computational and Graphical Statistics, 29(4):1–12.

P´olya, G. (1920). ¨Uber den zentralen grenzwertsatz der wahrscheinlichkeitsrechnung und das momentenproblem. Mathematische Zeitschrift, 8(3-4):171–181.

R Core Team (2018). R: A Language and Environment for Statistical Computing. R Foundation for Statistical Computing, Vienna, Austria.

Rivoirard, V. and Rousseau, J. (2012). Bernstein–von Mises theorem for linear functionals of the density. The Annals of Statistics, 40(3):1489–1523.

Rousseau, J. and Mengersen, K. (2011). Asymptotic behaviour of the posterior distribution in overfitted mixture models. Journal of the Royal Statistical Society: Series B (Statistical Methodology), 73(5):689–710.

Scott, D. W. (2015). Multivariate Density Estimation: Theory, Practice, and Visualization. John Wiley & Sons.

Sheather, S. J. and Jones, M. C. (1991). A reliable data-based bandwidth selection method for kernel density estimation. Journal of the Royal Statistical Society: Series B (Methodological), 53(3):683–690.

Song, H., Wang, Y., and Dunson, D. B. (2020). Distributed Bayesian clustering using finite mixture of mixtures. arXiv preprint, page arXiv:2003.13936.

Song, P. X.-K. (2000). Multivariate dispersion models generated from Gaussian copula. Scandinavian Journal of Statistics, 27(2):305–320.

Terrell, G. R. and Scott, D. W. (1992). Variable kernel density estimation. The Annals of Statistics, pages 1236–1265.

Teschl, G. (2009). Mathematical methods in quantum mechanics. Graduate Studies in Mathematics, 99:106.

Tsybakov, A. B. (2004). Introduction to nonparametric estimation, 2009. URL https://doi. org/10.1007/b13794. Revised and extended from the.

Vaidya, P. M. (1986). An optimal algorithm for the all-nearest-neighbors problem. In 27th Annual Symposium on Foundations of Computer Science, pages 117–122.

Wand, M. P. and Jones, M. C. (1994). Multivariate plug-in bandwidth selection. Computational Statistics, 9(2):97–116.

Wang, L. and Dunson, D. B. (2011). Fast Bayesian inference in Dirichlet process mixture models. Journal of Computational and Graphical Statistics, 20(1):196–216.

West, M. (1992). Hyperparameter estimation in Dirichlet process mixture models. Duke University ISDS Discussion Paper# 92-A03.

Wong, W. H. and Ma, L. (2010). Optional Polya tree and Bayesian inference. The Annals of Statistics, 38(3):1433–1459.

Zhang, X., Nott, D. J., Yau, C., and Jasra, A. (2014). A sequential algorithm for fast fitting of Dirichlet process mixture models. Journal of Computational and Graphical Statistics, 23(4):1143–1162.