Asymptotic analysis of the jittering kernel density estimator

# Asymptotic analysis of the jittering kernel density estimator

Thomas Nagler111Corresponding author, Department of Mathematics, Technische Universität München, Boltzmanstraße 3, DE-85748 Garching, Germany (email: thomas.nagler@tum.de)
July 6, 2019
###### Abstract

Abstract
Jittering estimators are nonparametric function estimators for mixed data. They extend arbitrary estimators from the continuous setting by adding random noise to discrete variables. We give an in-depth analysis of the jittering kernel density estimator, which reveals several appealing properties. The estimator is strongly consistent, asymptotically normal, and unbiased for discrete variables. It converges at minimax-optimal rates, which are established as a by-product of our analysis. To understand the effect of adding noise, we further study its asymptotic efficiency and finite sample bias in the univariate discrete case. Simulations show that the estimator is competitive on finite samples. The analysis suggests that similar properties can be expected for other jittering estimators.
Keywords: density, discrete, jittering, kernel, minimax, mixed data

## 1 Introduction

Multivariate density estimation is a central field in nonparametric statistics. Yet many popular methods have a significant drawback in applications: they can only be applied to continuous data. Some estimators have been specifically designed to allow for mixed continuous and discrete data (Ahmad and Cerrito, 1994, Li and Racine, 2003, Hall et al., 1983, Efromovich, 2011), but the number is small compared to the methods available in a purely continuous framework.

A common trick among practitioners is to make the discrete variables continuous by adding a small amount of noise. The noisy data is continuous and the usual nonparametric estimators apply. But the addition of random noise can introduce bias, so this procedure generally lacks justification. Nagler (2017) showed that adding noise still allows for valid estimates when the noise comes from a certain class of distributions. Then any nonparametric density estimator can be used in the mixed data setting. The resulting estimators are called jittering estimators.

Jittering estimators have so far been neglected in academic research, likely due to the widespread concern that jittering causes a loss in efficiency. The main objective of this article is to demonstrate that this concern is usually unjustified. To this end, we give an in-depth analysis of a simple instance from the class of jittering estimators: the jittering kernel density estimator, which is the jittering analog of the classical kernel density estimator (Parzen, 1962, Rosenblatt, 1956, Wand, 1992). We shall show that it maintains all the properties expected from a good nonparametric density estimator:

1. It is asymptotically normal and asymptotically unbiased for discrete variables (Theorem 1).

2. It is strongly and uniformly consistent (Theorem 2).

3. It is relatively efficient, even fully efficient in specific cases (Section 4.1).

4. It converges at minimax-optimal rates for a large class of target densities (Theorem 3 and Theorem 4). To the best of the author’s knowledge, these are the first results on minimax-optimality of nonparametric density estimators for mixed data.

Although focus is on only one instance of the class of jittering estimators, we can expect that others have similar properties.

The remainder of this article is organized as follows. Section 2 introduces the the jittering estimator and some assumptions. Section 3 gives a comprehensive asymptotic analysis which is complemented by a study of the asymptotic efficiency and finite sample bias in the univariate discrete setting (Section 4). Section 5 establishes with minimax-optimal rates for density estimation in a nonparametric mixed data model. Section 6 supports demonstrates that the estimator is also competitive on finite samples; Section 7 offers conclusions. Proofs of all theorems are deferred to Appendix A.

## 2 The estimator

Suppose that is a random vector with discrete component and continuous component . We explicitly allow for the cases where , (all variables are discrete) and , (all variables are continuous). Our goal is to estimate the density of based on ‘observations’ , , which are iid random vectors having the same distribution as . In this context, is the density with respect to the product of the counting and Lebesgue measures, i.e.,

 fZ,X(z,x)=∂q∂x1⋯∂xqPr(Z=z,X≤x).

Let be a real-valued function, called kernel, and abbreviate for any , . The classical kernel density estimator is defined as

 ˆf(z,x)=1nhpnbqnn∑i=1K(Zi−zhn)K(Xi−xbn), (1)

where are called bandwidth parameters and control the amount of smoothing. The above definition of the estimator is simplified to ease our exposition: we use only one parameter () for smoothing all components of and one parameter () for smoothing the components of . In practice, one would use a single parameter for each variable or even a bandwidth matrix (see, e.g., Scott, 2008).

The estimator only works for continuous random vectors. To make it applicable to mixed data, we make all discrete variables continuous by adding noise. Let , , be iid random vectors independent from , . Suppose further that the components of are iid with density . The jittering kernel density estimator is defined as the classical kernel density estimator applied to , :

 ˜f(z,x)=1nhpnbqnn∑i=1K(Zi+Ei−zhn)K(Xi−xbn). (2)

To facilitate our analysis, the following conditions are imposed on the kernel function:

###### Assumptions.

1. is a continuous function satisfying .

2. There is , , such that for ,

 ∫[0,1]tkK(t)dt=0,∫[0,1]tℓK(t)dt>0.
###### Remark 1.

A kernel function satisfying K2 is called -th order kernel (see, e.g., Marron, 1994). ∎

We further assume that the noise density belongs to the class , as defined in Nagler (2017):

###### Definition 1.

We say that for some , if

1. is an absolutely continuous probability density function,

2. for all ,

3. for all .

The density of is given by

 fη(z,x)=∑z′∈\mathdsZpf(z′,x)p∏j=1η(zj−z′j),(z,x)∈\mathdsZp×\mathdsRq.

The class ensures that is well-behaved. The most important properties are summarized in the following result (see, Nagler, 2017, Propositions 2 and 3).

###### Lemma 1.

Suppose the components of are iid with density . Then the joint density of satisfies for all , and such that ,

 fη(z,x)=f(z,x),∂¯¯¯¯mfη(z,x)∂zm11⋯∂zmpp=0.

The first equality implies that we can equivalently estimate instead of . This is convenient because is the density of a purely continuous random vector. The second equality states that all derivatives w.r.t.  vanish, which makes estimation even easier.

###### Remark 2.

The estimator is similar to the estimators of Ahmad and Cerrito 1994 and Li and Racine (2003). The difference lies in the kernel function for discrete data. The estimators of Ahmad and Cerrito 1994 and Li and Racine (2003) use a deterministic kernel function which is defined on the integers. In contrast, the jittering kernel density estimator (2) uses a random kernel defined on a compact subset of where randomness is induced by .

## 3 Asymptotic analysis in the general setting

### 3.1 Asymptotic distribution

We first study the asymptotic distribution of the jittering kernel density estimator. To motivate our first theorem, we recall a a classical result from kernel density estimation in the purely continuous setting (e.g., Wand, 1992). If is the density of a continuous random vector , sufficiently smooth, , and , , then

 (3)

Recall that is nothing else than applied to , . Lemma 1 showed that , the density of , has vanishing derivatives with respect to . We can thus expect the first sum in the bias term in (3) to vanish asymptotically. In fact, it becomes exactly zero when . The following result improves upon the properties implied by (3) by taking these considerations into account.

###### Assumptions.

1. is times continuously differentiable with respect to .

2. and hold with .

3. .

4. and as .

5. There is , such that for all .

###### Theorem 1.

Under assumptions A1-A5, it holds for any ,

 E{˜f(z,x)} =f(z,x)+bℓnσℓℓ!q∑j=1∂ℓf(z,x)∂xℓj+o(bℓn), Var{˜f(z,x)} =f(z,x)nbqn{h−pnκp+q−bqnf(z,x)}+o(1nhpnbqn),

where and . If further ,

 ˜f(z,x)−E{˜f(z,x)}Var{˜f(z,x)}\lx@stackreld→N(0,1).
###### Remark 3.

The assumptions in Theorem 1 differ from those usually made in the continuous framework. There are no assumptions on the smoothness of with respect to , because its local behavior is controlled by . Further, is not required to vanish asymptotically, but should be less than for large . This is sufficient to ensure that there is no bias with respect to . Further decreasing does not change the bias, but inflates the variance. ∎

###### Remark 4.

The asymptotic variance does not involve on or its class parameters and (and neither does the asymptotic bias). Intuitively, we would expect an increase in the estimator’s variance because we are adding random noise. Apparently this effect is dominated by the sampling variability in the original data and asymptotically negligible. So there should be no benefit from averaging over multiple jitters (at least asymptotically). This is in contrast to empirical processes of jittered data (Genest et al., 2017).

### 3.2 Asymptotically optimal bandwidths

A standard tool for studying optimal bandwidths is the asymptotic mean squared error,

Under the assumptions of Theorem 1, we get

 AMSE{˜f(z,x)}

For , it is easy to check that the bandwidth minimizing the AMSE satisfies . This is well-known as the optimal rate for the classical kernel density estimator when . The AMSE further suggests that it is optimal to choose as large as possible. The largest allowed by A5 is . Asymptotically, this is the optimal bandwidth. We shall see shortly that this choice means that we are not smoothing the discrete variables at all. This is not unreasonable: in contrast to the continuous case, smoothing discrete variables is not necessary for consistent nonparametric estimation (for a discussion, see, Simar et al., 2011).

On finite samples can be too small. If , the estimator can be written as

 ˜f(z,x)=1nhpnbqn∑i:Zi=zK(Eihn)K(Xi−xbn),

Indeed, the estimator neglects all observations where and, thus, does not smooth with respect to the discrete variables. This also means that if for all . Theorem 1 implicitly assumes that is large enough to provide sufficiently many observations with . This is guaranteed asymptotically whenever ), but often demands sample sizes much larger than what is common.

We conclude that Theorem 1 is not useful for bandwidth selection on samples of small or moderate size. Cross-validation techniques are more appropriate tools in the mixed data setting (see, e.g., J. Aitchison, 1976, Racine and Li, 2004, Hall et al., 2004).

### 3.3 Consistency

Theorem 1 implies pointwise consistency of the jittering kernel density estimator, but assumption A1 is more strict than necessary. The following result weakens this assumption and additionally establishes strong uniform consistency.

###### Assumptions.

1. The th derivative of exists and is uniformly Lipschitz on .

###### Theorem 2.

Suppose that assumptions , hold. Then, for all ,

 ˜f(z,x)−f(z,x) =Op{bℓn+(nhpnbqn)−1/2}, (4) supS∣∣˜f(z,x)−f(z,x)∣∣ =Oa.s.{bℓn+(max{lnlnn,lnh−1n,lnb−1n}nhpnbqn)1/2}. (5)
###### Remark 5.

If there are such that for all , the rates of convergence in Theorem 2 do not involve , the dimension of the discrete variables. So adding more discrete variables does not change the convergence rate of the estimator. In particular, there is no cost for recoding unordered categorical variables into several binary variables. ∎

###### Remark 6.
1. The best rate in (4) is and achieved when and .

2. For , the best rate in (5) is and achieved when , .

3. For , the best rate in (5) is and achieved when .

## 4 A closer look at the univariate discrete setting

The jittering kernel density estimator handles continuous variables just like the classical kernel density estimator. How it smooths discrete variables is less obvious. To gain a better understanding, we study its asymptotic efficiency and finite sample bias when there is only one discrete variable (, ).

### 4.1 Asymptotic efficiency

For convenience, set . The expectation and variance in Theorem 1 become

 E{˜f(z)}=f(z),Var{˜f(z)}=f(z)n[min(γ1,1−γ2)−1κ−f(z)]+o(n−1),

The most efficient point estimator for a discrete probability is the sample frequency . It satisfies

 E{fn(z)}=f(z),Var{fn(z)}=f(z)n{1−f(z)}.

The asymptotic relative efficiency (ARE) of relative to is defined as

 ARE{˜f(z):fn(z)}=AVar{fn(z)}AVar{˜f(z)},

where denotes the leading term of an asymptotic expansion of the variance. The ARE is interpreted as follows: If the estimator is used with observations, then one needs observations to obtain the same accuracy with . If the ARE is less than one, then needs less observations, i.e., is more efficient than . If the ARE is greater then one, it is the other way around. If it is exactly one, the two estimators are equally efficient.

Straightforward calculations yield

 ARE{˜f(z):fn(z)} =1−f(z)min{γ1,1−γ2}−1κ−f(z) =(1+min{γ1,1−γ2}−1κ−11−f(z))−1≤1.

The relative efficiency depends on three quantities:

• It is increasing in and the most efficient choice is , which corresponds to the uniform error density on . On the other hand, the relative efficiency approaches 0 for or .

• It is decreasing in , which is the roughness of the kernel . The ‘least rough‘ kernel is the is the uniform kernel, i.e., , for which . But this kernel is rather unpopular in practice. A more widely used kernel is the Epanechnikov kernel, , for which .

• It is decreasing in . The worst case is that , for which the ARE is zero. For a variable, noise, and the Epanechnikov kernel, we get .

###### Remark 7.

Suppose is the uniform density on (for which ), , and is the uniform kernel (for which ). Then, the two estimators are equally efficient. In fact, the estimator becomes

 ˜f(z) =1nhnn∑i=12−1\mathds1(|Zi+Ei−z|≤hn) =2nn∑i=12−1\mathds1(|Zi+Ei−z|≤1/2) =1nn∑i=1\mathds1(Zi=z),

which is exactly the sample frequency estimator . ∎

### 4.2 Finite sample bias

Assuming , Theorem 1 shows that is unbiased in a purely discrete setting. On small samples, it is often necessary to choose a larger bandwidth (see Section 3.2). When , the estimator is usually biased.

###### Lemma 2.

Suppose that and satisfies . Then,

 E{˜f(z)}−f(z) = ⌈hn−1/2⌉∑k=1ρηk(hn)f(z+k)−{ρηk(hn)+ρη−k(hn)}f(z)+ρη−k(hn)f(z−k)k2,

where and

 Aηk(hn)=[(1−γ2−k)h−1n,(−1+γ2−k)h−1n]∩[−1,1].

To interpret the bias, it is helpful to focus on a simple case first.

###### Corollary 1.

Suppose that and is a symmetric function satisfying . Then for all ,

 E{˜f(z)}=f(z)+⌈hn−1/2⌉∑k=1ρk(hn)Δ2kf(z),

where

 Δ2kf(z)=f(z+k)−2f(z)+f(z−k)k2,

and with

The operator is known as the second order central difference operator (e.g., Monahan, 2011). It is commonly used as numerical approximation of second order derivative of real-valued functions, which is

 d2f(x)dx2=lims→0f(x+s)−2f(x)+f(x−s)s2.

We can interpret as a discrete analogue to the second order derivative of a real-valued function. In this aspect, the discrete setting is similar to the continuous one (where the bias of is proportional to the second order derivative).

The parameter is called the step size and determines how local the derivative approximation is. The bias of is a weighted sum of such ‘derivatives’ for several values of . The bandwidth limits the maximal step size and thereby controls the locality of the bias. Although not universally true, smaller values of typically correspond to a smaller bias. A simple counter example is when for all , where the bias is zero for all . There are also situations where decreasing leads to a larger bias. This phenomenon also exists in the continuous setting, but is disguised by asymptotic approximations. When as in Theorem 1, the estimator is unbiased.

The bias in Lemma 2 can be interpreted similarly. But is replaced by a weighted approximation of the derivative. If or are asymmetric, different weights will be assigned to the ‘forward derivative’ and the ‘backward derivative’ .

## 5 Minimax rate optimality

The maximum risk associated with a class of densities and a (semi-) distance is defined as

 Rn(ˆf,F,d)=supf∈FEf{d2(ˆf,f)}, (6)

We consider two semi-distances that relate to pointwise and uniform consistency of , respectively:

 d(z,x)(ˆf,f)= ∣∣ˆf(z,x)−f(z,x)∣∣,for some (z,x)∈\mathdsZp×\mathdsRq, d∞,S(ˆf,f)= supS∣∣ˆf(z,x)−f(z,x)∣∣,for some S⊂\mathdsZp×\mathdsRq.

For , we shall consider all bounded density functions whose continuous part belongs to a Hölder class. For , we use the multi-index notations , and denote the partial derivatives of with respect to as

 Daxf(z,x)=∂|a|f(z,x)∂a1x1⋯∂aqxq. (7)
###### Definition 2.

For and , , , the class is defined as all functions such that for all with ,

1. is a probability density on ,

2. exists for all and

 supz∈\mathdsZp,x,x′∈\mathdsRq{∣∣Daxf(z,x)−Daxf(z,x′)∣∣∥x−x′∥α2+f(z,x)}≤λ
###### Remark 8.

If and , contains all densities on . If and , it is a Hölder class on . ∎

The following result establishes convergence rates of the jittering kernel density estimator with respect to the maximum risk.

###### Theorem 3.

Denote as the estimator defined in Equation 2. Suppose and assumptions A2–A4 of Theorem 1 hold with , , , . Assume further that there are such that for all . Then there exists such that

 limsupn→∞r−2nRn(ˆf,F,d)≤¯¯c,

in each of the following cases:

1. ,

2. , , ,

3. , , , ,

4. , , , ,

for arbitrary and .

We shall see that the rates in Theorem 3 (i)–(iii) are optimal in a minimax sense. The minimax risk is defined as

where the infimum is taken over all possible estimators of . In our context, an ‘estimator’ is any measurable function of , .

###### Definition 3.

A sequence of positive real numbers is called

1. an upper bound on the minimax rate if there is such that

 limsupn→∞r−2nR∗n(F,d)≤¯¯c.
2. a lower bound on the minimax rate if there is such that

 liminfn→∞r−2nR∗n(F,d)≥c–,
3. a minimax-optimal rate of convergence if both (i) and (ii) hold.

In a purely continuous setting, optimal rates have long been established (Stone, 1980, 1983, Ibragimov and Khas’ minskii, 1983). To the best of the author’s knowledge, there are no results on optimal rates in the mixed data setting.

To show that a rate is minimax-optimal, we have to check that it is both an upper and lower bound on the minimax rate. Theorem 3 already gives us an upper bound, since, for any estimator ,

 R∗n(F,d)=infˆfRn(ˆf,F,d)≤Rn(ˆf,F,d).

Lower bounds on the minimax rate can be deduced easily by considering subsets of for which lower bounds are known (see Section A.4).

###### Theorem 4.

Let and . The minimax-optimal rate of convergence associated with the class and distance satisfies

1. , for

2. , for , ,

3. , for , , ,

4. , for , , ,

###### Remark 9.

Theorem 3 and Theorem 4 imply that the jittering kernel density estimator converges at minimax-optimal rates for cases (i)–(iii). ∎

###### Remark 10.

Theorem 4 only provides an interval for the optimal rate in case (iv). Minimax analysis for this setting is surprisingly har; see (Han et al., 2015) for minimax rates with respect to the distance. The interval is quite narrow, differing only by a factor of size . The exact rate, however, remains an open problem. ∎

## 6 Simulation experiments

The jittering kernel density estimator has appealing asymptotic properties. This may come as a surprise: since we are adding noise to the data, we could expect that the data become less informative and uncertainty increases. We complement our asymptotic arguments with a small numerical experiment that illustrates the small sample performance of the estimator. Because of its wide use and close resemblance to our approach, we will use the estimator of Li and Racine (2003) as a benchmark.

We use the following setup:

• We compare three estimators

1. jkde: the jittering kernel density estimator with noise density , for which .

2. jkde2: the jittering kernel density estimator with noise density (as in, Nagler, 2017, Example 3), for which , .

3. liracine: the estimator of Li and Racine (2003) as implemented in the np package (Hayfield and Racine, 2008).

Contrary to (2), we use one bandwidth parameter for each variable. Both estimators use likelihood cross-validation for bandwidth selection.

• We estimate the density of a vector , where for all , for all . For sake of simplicity, all variables are simulated independently.

• Results are based on simulated data sets with sample sizes .

• As a performance measure we use the root average square error (RASE) computed over a grid in . More specifically, we use , , and

 RASE(ˆf,f)=√∑z1∈Z⋯∑zp∈Z∑x1∈X⋯∑xq∈X{ˆf(z,x)−f(z,x)}2. Figure 1: RASE achieved by the two estimators for various choices of p, q, and m. Each estimator is represented by two boxes; the left box corresponds to n=50, the right to n=200.

Figure 1 shows the estimators’ performance for various values of , and . Each estimator is represented by two boxes, where the left box corresponds to and the right box to . The choice of noise density seems to be of minor importance: jkde and jkde2 give almost identical results. Compared to liracine, the two estimator show only subtle differences. The two jittering estimators are more accurate in all scenarios with , and less accurate when . This is related to our observation from Section 4.1 that the efficiency is worse when is large. The relative performance of the three estimators is consistent across the two sample sizes under consideration. Overall, the jittering estimators are competitive with the benchmark estimator liracine. We found no evidence that adding artificial noise negatively affects the accuracy of the estimates. This confirms what was suggested by the estimator’s asymptotic properties.

## 7 Conclusion

This article gave an in-depth analysis of the behavior of the jittering kernel density estimator. It was shown to have appealing large-sample properties and perform well on small samples.

Although our focus was on a particular instance of the class of jittering estimators, we also learned something about the class as a whole. Adding noise to discrete variables does not have a negative impact on estimation accuracy. This is true for both large samples (as confirmed by our asymptotic analysis) and small samples (as illustrated by simulations). More specifically, it allows for estimators that are optimal in terms of convergence rates and efficiency. It is likely that these findings generalize to more sophisticated density estimators or estimators of functionals of the density, such as regression functions.

### Acknowledgements

This work was partially supported by the German Research Foundation (DFG grant CZ 86/5-1).

## Appendix A Proofs

### a.1 Proof of Theorem 1

We first calculate the bias term. Using a change of variables, we get

 E{˜f(z,x)} =1hpnbqnE{K(Z+E−zhn)K(X−xbn)} =1hpnbqn∫\mathdsRp+qK(s−zhn)K(t−xbn)fη(s,t)dsdt =∫\mathdsRp+qK(u)K(v)fη(z+hnu,x+bnv)dudv

Since , it holds for all and that . Furthermore, is zero outside of . Hence, for ,

 E{˜f(z,x)} =∫\mathdsRp+qK(u)K(v)f(z,x+bnv)dudv =∫[−1,1]qK(v)f(z,x+bnv)dv. (8)

Recall the derivative notation from (7). An -th order Taylor expansion of yields that

for some , where the second equality is due to K2. The second sum is because all terms are bounded by A1 and K1. In summary,

 E{˜f(z,x)}−f(z,x) =bℓnσℓℓ!q∑j=1∂ℓf(z,x)∂xℓj+o(bℓn),

as claimed.

For the variance, we get

 Var{˜f(z,x)} =1n[1h2pnb2qnE{K(Z+E−zhn)2K(X−xbn)2} =1n[−1h2pnb2qnE{K(Z+E−zhn)K(X−xbn)}2].

The second term in square brackets has already been calculated for the bias. Using similar arguments, we can show

 =1nhpnbqn∫\mathdsRp+qK2(s−zhn)K2(t−xbn)fη(s,t)dsdt =κp∫[−1,1]qK2(v)f(z,x+bnv)dv =κp+qf(z,x)+o(1).

Together,