Ake: An R Package for Discrete and Continuous Associated Kernel Estimations

1 Introduction

Kernel smoothing methods are popular tools for revealing the structure of data that could be missed by parametric methods. For real datasets, we often encounter continuous (bounded, positive) or discrete (count, categorical) data types. The classical kernels methods assume that the underlying distribution is unbounded continuous, which is frequently not the case; see, for example, (Duong 2007) for multivariate kernel density estimation and discriminant analysis. A solution is provided for categorical data sets by (Hayfield and Racine 2008). In fact, they used kernels well adapted for these categorical sets (Aitchison and Aitken 1976). Throughout the present paper, the unidimensional support $\mathbb{T}$ of the variable of interest can be $\{0,1,\dots ,N\}$, $[a,b]$ or $[0,\mathrm{\infty })$ for a given integer $N$ and reals $a<b$.

The recently developed Ake package, implements associated kernels that seamlessly deal with continuous (bounded, positive) and discrete (categorical, count) data types often found in applied settings; see, for example, (Libengué 2013) and (Kokonendji and Senga Kiessé 2011). These associated kernels are used to smooth probability density functions (p.d.f.), probability mass functions (p.m.f.) or regression functions. The coming versions of this package will contain, among others, p.d.f. estimation of heavy tailed data (e.g., Ziane et al. 2015) and the estimation of other functionals. The bandwidth selection remains crucial in associated kernel estimations of p.d.f., p.m.f. or regression functions. Some methods have been investigated for selecting bandwidth parameters but the commonly used is the least squared cross-validation. A Bayesian approach has been also recently introduced by (Zougab et al. 2012) in the case of a binomial kernel. This method can be extended to various associated kernels with other functionals. Despite the great number of packages implemented for nonparametric estimation in continuous cases with unbounded kernels, to the best of our knowledge, the R packages to estimate p.m.f. with categorical or count variables, p.d.f. with bounded or positive datasets, and regression functions have been far less investigated.

The rest of the paper is organized as follows. In Section 2, we briefly describe the definition of associated kernels and then illustrate examples in both continuous and discrete cases which are discussed. Then, the associated kernel estimator for p.d.f. or p.m.f. is presented and illustrated with some R codes in Section 3. In particular, three bandwidth selection methods are available: cross-validation for any (continuous or discrete) associated kernel, the Bayesian local method for the binomial kernel and also a new theoretical Bayesian adaptive method for the gamma kernel. Also, some practical and theoretical aspects of the normalizing constant in both p.d.f. and p.m.f. estimations are given. Section 4 investigates the case of regression functions with two bandwidth selection techniques: cross-validation and also the Bayesian global method for the binomial kernel. Section 5 concludes.

2 Non-classical associated kernels

Recall that the support $\mathbb{T}$ of the p.m.f., p.d.f. or regression function, to be estimated, is any set $\{0,1,\dots ,N\}$, $[a,b]$ or $[0,\mathrm{\infty })$ for a given integer $N$ and reals $a<b$. The associated kernel in both continuous and discrete cases is defined as follows.

Definition 2.2.1

(Kokonendji and Senga Kiessé 2011),(Libengué 2013) Let $\mathbb{T}(\subseteq \mathbb{R})$ be the support of the p.m.f., p.d.f. or regression function, to be estimated, $x\in \mathbb{T}$ a target and $h$ a bandwidth. A parametrized p.m.f. (respectively p.d.f.) $K_{x,h}(\cdot )$ of support ${\mathbb{S}}_{x,h}(\subseteq \mathbb{R})$ is called “associated kernel” if the following conditions are satisfied:

$$\begin{array}{}\text{(1)}& \begin{array}{rl}& x\in {\mathbb{S}}_{x,h},\end{array}\end{array}$$

$$\begin{array}{}\text{(2)}& \begin{array}{rl}& \mathbb{E}\left({\mathcal{Z}}_{x,h}\right)=x+a(x,h),\end{array}\end{array}$$

$$\begin{array}{}\text{(3)}& \begin{array}{rl}& \mathbb{V}\left({\mathcal{Z}}_{x,h}\right)=b(x,h),\end{array}\end{array}$$

where ${\mathcal{Z}}_{x,h}$ denotes the random variable with p.m.f. (respectively p.d.f.) $K_{x,h}$ and both $a(x,h)$ and $b(x,h)$ tend to $0$ as $h$ goes to $0$.

Remark 2.2.2

This definition has the following interesting interpretations:

1. The function $K_{x,h}(\cdot )$ is not necessary symmetric and is intrinsically linked to $x$ and $h$.

2. The support ${\mathbb{S}}_{x,h}$ is not necessary symmetric around $x$; it can depend or not on $x$ and $h$.

3. The condition ((1)) can be viewed as ${\cup }_{x\in \mathbb{T}}{\mathbb{S}}_{x,h}\supseteq \mathbb{T}$ and it implies that the associated kernel takes into account the support $\mathbb{T}$ of the density $f$, to be estimated.

4. If ${\cup }_{x\in \mathbb{T}}{\mathbb{S}}_{x,h}$ does not contain $\mathbb{T}$ then this is the well-known problem of boundary bias.

5. Both conditions ((2)) and ((3)) indicate that the associated kernel is more and more concentrated around $x$ as $h$ goes to $0$. This highlights the peculiarity of the associated kernel which can change its shape according to the target position.

In order to construct an associated kernel $K_{x,h}(\cdot )$ from a parametric (discrete or continuous) probability distribution $K_{\theta }$, $\theta \in \mathrm{\Theta }\subset {\mathbb{R}}^d$ on the support ${\mathbb{S}}_{\theta }$ such that ${\mathbb{S}}_{\theta }\cap \mathbb{T}\ne \mathrm{\varnothing }$, we need to establish a correspondence between $(x,h)\in \mathbb{T}\times (0,\mathrm{\infty })$ and $\theta \in \mathrm{\Theta }$; see (Kokonendji and Senga Kiessé 2011). In what follows, we will call $K\equiv K_{\theta }$ the type of kernel to make a difference from the classical notion of a continuous symmetric (e.g., Gaussian) kernel. In this context, the choice of the associated kernel becomes important as well as that of the bandwidth. Moreover, we distinguish the associated kernels said sometimes of “second order” of those said of “first order” which verify the two first conditions (1) and (2). The rest of this section is devoted to discuss examples of associated kernels in both discrete and continuous cases.

Discrete associated kernels

Among the discrete associated kernels found in literature, we here use the best in sense of Definition 2.2.1. Negative binomial and Poisson kernels are respectively overdispersed (i.e., $\mathbb{V}({\mathcal{Z}}_{x,h})>\mathbb{E}({\mathcal{Z}}_{x,h}))$ and equisdispersed (i.e., $\mathbb{V}({\mathcal{Z}}_{x,h})=\mathbb{E}({\mathcal{Z}}_{x,h}))$ and thus are not recommended; see (Kokonendji and Senga Kiessé 2011) for further details. The first associated kernel listed below, namely the binomial kernel, is the best of the first order or standard kernels which satisfies $$\begin{array}{}\text{(4)}& \underset{h\to 0}{lim}\mathbb{V}({\mathcal{Z}}_{x,h})\in \mathcal{V}(0),\end{array}$$ where $\mathcal{V}(0)$ is a neighborhood of $0$ which does not depend on $x$. The two other discrete associated kernels satisfy all conditions of Definition 2.2.1.

• The binomial (bino) kernel is defined on the support ${\mathbb{S}}_x=\{0,1,\dots ,x+1\}$ with $x\in \mathbb{T}:=\mathbb{N}=\{0,1,\dots \}$ and then $h\in (0,1]$: $$B_{x,h}(u)=\frac{(x+1)!}{u!(x+1-u)!}{\left(\frac{x+h}{x+1}\right)}^u{\left(\frac{1-h}{x+1}\right)}^{x+1-u}{\mathbb{1}}_{{\mathbb{S}}_x}(u),$$ where ${\mathbb{1}}_A$ denotes the indicator function of any given event $A$. Note that $B_{x,h}$ is the p.m.f. of the binomial distribution $\mathcal{B}(x+1;(x+h)/(x+1))$ with its number of trials $x+1$ and its success probability in each trial $(x+h)/(x+1)$. It is appropriate for count data with small or moderate sample sizes and, also, it satisfies (4) rather than ((3)); see  (Kokonendji and Senga Kiessé 2011) and also (Zougab et al. 2012) for a bandwidth selection by Bayesian method.

• The following class of symmetric discrete triangular kernels has been proposed in (Kokonendji et al. 2007). The support $\mathbb{T}$ of the p.m.f. $f$ to be estimated, can be unbounded (e.g., $\mathbb{N},\mathbb{Z}$) or finite (e.g., $\{0,1,\dots ,N\}$). Then, suppose that $h$ is a given bandwidth parameter and $a$ is an arbitrary and fixed integer. For fixed arm $a\in \mathbb{N}$, the discrete triangular (DTra) kernel is defined on ${\mathbb{S}}_{x,a}=\{x,x\pm 1,\dots ,x\pm a\}$ with $x\in \mathbb{T}=\mathbb{N}$: $$D{T}_{x,h;a}(u)=\frac{(a+1{)}^h-|u-x{|}^h}{P(a,h)}{\mathbb{1}}_{{\mathbb{S}}_{x,a}}(u),$$ where $P(a,h)=(2a+1)(a+1)-2\sum _{k=0}^a{k}^h$ is the normalizing constant. It is symmetric around the target $x$, satisfying Definition 2.2.1 and suitable for count variables; see  (Kokonendji and Zocchi 2010) for an asymmetric version. Note that $h\to 0$ gives the Dirac kernel.

• A discrete kernel estimator for categorical data has been introduced in (Aitchison and Aitken 1976). Its asymmetric discrete associated kernel version that we here label DiracDU (DirDU) as “Dirac Discrete Uniform” has been deduced in (Kokonendji and Senga Kiessé 2011) as follows. For fixed $c\in \{2,3,\dots \}$ the number of categories, we define ${\mathbb{S}}_c=\{0,1,\dots ,c-1\}$ and $$D{U}_{x,h;c}(u)=(1-h){\mathbb{1}}_{\left\{x\right\}}(u)+{\displaystyle \frac{h}{c-1}}{\mathbb{1}}_{{\mathbb{S}}_c\setminus \left\{x\right\}}(u),$$ where $h\in (0,1]$ and $x\in \mathbb{T}={\mathbb{S}}_c$. In addition, the target $x$ can be considered as the reference point of $f$ to be estimated; and, the smoothing parameter $h$ is such that $1-h$ is the success probability of the reference point. This DiracDU kernel is symmetric around the target, satisfying Definition 2.2.1 and appropriated for categorical set $\mathbb{T}$. See, e.g., (Racine and Li 2004) for some uses. Note that $h=0$ provides the Dirac kernel.

Continuous associated kernels

One can find several continuous associated kernels in literature among the Birnbaum-Saunders of (Jin and Kawczak 2003). Here, we present seven associated kernels well adapted for the estimations of density or regression functions on any compact or nonnegative support of datasets. All these associated kernels satisfy Definition 2.2.1.

• The extended beta (BE) kernel is defined on ${\mathbb{S}}_{x,h,a,b}=[a,b]=\mathbb{T}$ with $a<b<\mathrm{\infty }$, $x\in \mathbb{T}$ and $h>0$ such that

  $$B{E}_{x,h,a,b}(u)=\frac{(u-a{)}^{(x-a)/\{(b-a)h\}}(b-u{)}^{(b-x)/\{(b-a)h\}}}{(b-a{)}^{1+h^{-1}}B(1+(x-a)/(b-a)h,1+(b-x)/(b-a)h)}{\mathbb{1}}_{{\mathbb{S}}_{x,h,a,b}}(u),$$ where $B(r,s)={\int }_0^1{t}^{r-1}(1-t{)}^{s-1}dt$ is the usual beta function with $r>0$, $s>0$; see (Libengué 2013). For $a=0$ and $b=1$, it corresponds to the beta kernel (Chen 1999) which is the p.d.f. of the beta distribution with shape parameters $1+x/h$ and $(1-x)/h$. The extended beta kernel is appropriate for any compact support of observations.

• The gamma (GA) kernel is given on ${\mathbb{S}}_{x,h}=[0,\mathrm{\infty })=\mathbb{T}$ with $x\in \mathbb{T}$ and $h>0$: $$\begin{array}{}\text{(5)}& G{A}_{x,h}(u)={\displaystyle \frac{u^{x/h}}{\mathrm{\Gamma }(1+x/h)h^{1+x/h}}}\exp (-{\displaystyle \frac{u}{h}}){\mathbb{1}}_{[0,\mathrm{\infty })}(u),\end{array}$$ where $\mathrm{\Gamma }(v)={\int }_0^{\mathrm{\infty }}{s}^{v-1}\exp (-s)ds$ is the classical gamma function with $v>0$; see (Chen 2000). It is the p.d.f. of the gamma distribution $\mathcal{GA}(1+x/h,h)$ with scale parameter $1+x/h$ and shape parameter $h$. It is suitable for the non-negative real set $\mathbb{T}=[0,\mathrm{\infty })$.

• The lognormal (LN) kernel is defined on ${\mathbb{S}}_{x,h}=[0,\mathrm{\infty })=\mathbb{T}$ with $x\in \mathbb{T}$ and $h>0$ such that

  $$L{N}_{x,h}(u)=\frac{1}{uh\sqrt{2\pi }}\exp \{-\frac{1}{2}{(\frac{1}{h}\log (\frac{u}{x})-h)}^2\}{\mathbb{1}}_{{\mathbb{S}}_{x,h}}(u);$$ see (Libengué 2013) and also (Igarashi and Kakizawa 2015). It is the p.d.f. of the classical lognormal distribution with mean $\log (x)+h^2$ and standard deviation $h$.

• The reciprocal inverse Gaussian (RIG) kernel is given on ${\mathbb{S}}_{x,h}=(0,\mathrm{\infty })=\mathbb{T}$ with $x\in \mathbb{T}$ and $h>0$:

  $$RI{G}_{x,h}(u)=\frac{1}{\sqrt{2\pi hu}}\exp \{-\frac{(x^2+xh{)}^{1/2}}{2h}(\frac{u}{(x^2+xh{)}^{1/2}}-2+\frac{(x^2+xh{)}^{1/2}}{u})\}{\mathbb{1}}_{{\mathbb{S}}_{x,h}}(u);$$

  see (Scaillet 2004), (Libengué 2013) and also (Igarashi and Kakizawa 2015). It is the p.d.f. of the classical reciprocal inverse Gaussian distribution with mean $1/\sqrt{x^2+xh}$ and standard deviation $1/h$.

Remark 2.2.3

The three continuous associated kernels inverse gamma, inverse Gaussian and Gaussian are not adapted for density estimation on supports $[0,\mathrm{\infty })$ and thus are not included in the Ake package; see Part (b) of Figure 1.

Indeed:

• The inverse gamma (IGA) kernel, defined on ${\mathbb{S}}_{x,h}=(0,\mathrm{\infty })=\mathbb{T}$ with $x\in (0,1/h)$ and $h>0$ such that $$IG{A}_{x,h}(u)={\displaystyle \frac{h^{1-1/(xh)}}{\mathrm{\Gamma }(-1+1/(xh))}}{u}^{-1/(xh)}\exp (-{\displaystyle \frac{1}{hu}}){\mathbb{1}}_{(0,\mathrm{\infty })}(u)$$ (Libengué 2013), is graphically the worst since it does not well concentrate on the target $x$. Note that it is the p.d.f. of the inverse gamma distribution with scale parameter $-1+1/(xh)$ and scale parameter $1/h$.

• Also, the inverse Gaussian (IG) kernel, defined on ${\mathbb{S}}_{x,h}=(0,\mathrm{\infty })=\mathbb{T}$ with $x\in (0,1/3h)$ and $h>0$ by $$I{G}_{x,h}(u)=\frac{1}{\sqrt{2\pi hu}}\exp \{-\frac{(1-3xh{)}^{1/2}}{2h}(\frac{u}{(1-3xh{)}^{1/2}}-2+\frac{(1-3xh{)}^{1/2}}{u})\}{\mathbb{1}}_{{\mathbb{S}}_{x,h}}(u)$$ (Scaillet 2004),(Libengué 2013), has the same graphical properties as the inverse gamma. Note that it is the p.d.f. of the inverse Gaussian distribution $\mathcal{IG}(1+x/h,h)$ with scale parameter $x/(1-3xh{)}^{1/2}$ and shape parameter $1/h$.

• From the well known Gaussian kernel $K^G(u)=(h\sqrt{2\pi }{)}^{-1}\exp (u^2){\mathbb{1}}_{\mathbb{R}}(u)$, we define its associated version (Gaussian) on ${\mathbb{S}}_{x,h}=\mathbb{R}$ with $x\in \mathbb{T}:=\mathbb{R}$ and $h>0$: $$K_{x,h}^G(u)=\frac{1}{h\sqrt{2\pi }}\exp \left\{\frac{1}{2}{\left(\frac{u-x}{h}\right)}^2\right\}{\mathbb{1}}_{\mathbb{R}}(u).$$ It has the same shape at any target and thus is well adapted for continuous variables with unbounded supports but not for $[0,\mathrm{\infty })$ or compact set of $\mathbb{R}$; see also (Epanechnikov 1969) for another example of a continuous symmetric kernel.

Figure 1: Shapes of univariate (discrete and continuous) associated kernels: (a) DiracDU, discrete triangular a = 3 and binomial with same target x = 3 and bandwidth h = 0.13; (b) lognormal, inverse Gaussian, gamma, reciprocal inverse Gaussian, inverse gamma and Gaussian with same target x = 1.3 and h = 0.2.

Figure 1 shows some forms of the above-mentioned univariate associated kernels. The plots highlight the importance given to the target point and around it in discrete (a) and continuous (b) cases. Furthermore, for a fixed bandwidth $h$, the Gaussian keeps its same shape along the support; however, they change according to the target for the other non-classical associated kernels. This explains the inappropriateness of the Gaussian kernel for density or regression estimation in any bounded interval and of the DiracDU kernel for count regression estimation; see Part (b) of Figure 1. From Part (v) of Remark 2.2.2, the inverse gamma and inverse Gaussian are the worst since they do not well concentrate on the target $x$; see Remark 2.2.3. These previous associated kernels can be applied to various functionals.

We have implemented in R the method kern.fun for both discrete and continuous associated kernels. Seven possibilities are allowed for the kernel function. We enumerate the arguments and results of the default kern.fun.default function in Table 1. The kern.fun is used as follows for the binomial kernel:

R> x <- 5
R> h <- 0.1
R> y <- 0:10
R> k_b <- kern.fun(x, y, h, "discrete", "bino")

Table 1: Summary of arguments and results of kern.fun.

Arguments	Description
x	The target.
t	The single or the grid value where the function is computed.
h	The bandwidth or smoothing parameter.
ker	The associated kernel.
a0,a1	The left and right bounds of the support for the extended beta kernel.
a	The arm of the discrete triangular kernel. Default value is 1.
c	The number of categories in DiracDU kernel. Default value is 2.
Result	Description
	Returns a single value of the associated kernel function.

3 Density or probability mass function estimations

The p.d.f. or p.m.f. estimation is an usual application of the associated kernels. Let $X_1,\dots ,X_n$ be independent and identically distributed (i.i.d.) random variables with an unknown p.d.f. (respectively p.m.f.) $f$ on $\mathbb{T}$. An associated kernel estimator ${\hat{f}}_n$ of $f$ is simply: $$\begin{array}{}\text{(6)}& {\hat{f}}_n(x)={\displaystyle \frac{1}{n}}{\displaystyle \sum _{i=1}^n{K}_{x,h}(X_i),x\in \mathbb{T}.}\end{array}$$ Here, we point out some pointwise properties of the estimator (6) in both discrete and continuous cases.

Proposition 2.3.1

(Kokonendji and Senga Kiessé 2011),(Libengué 2013) Let $X_1,X_2,\dots ,X_n$ be an $n$ random sample i.i.d. from the unknown p.m.f. (respectively p.d.f.) $f$ on $\mathbb{T}$. Let ${\hat{f}}_n={\hat{f}}_{n,h,K}$ be an estimator (6) of $f$ with an associated kernel. Then, for all $x\in \mathbb{T}$ and $h>0$, we have $$\mathbb{E}\{\hat{f}(x)\}=\mathbb{E}\{f({\mathcal{Z}}_{x,h})\},$$ where ${\mathcal{Z}}_{x,h}$ is the random variable associated to the p.m.f. (respectively p.d.f.) $K_{x,h}$ on ${\mathbb{S}}_{x,h}$. Furthermore, for a p.m.f. (respectively p.d.f.), we have respectively ${\hat{f}}_n(x)\in [0,1]$ (respectively ${\hat{f}}_n(x)>0$) for all $x\in \mathbb{T}$ and $$\begin{array}{}\text{(7)}& {\int }_{x\in \mathbb{T}}{\hat{f}}_n(x)\nu (dx)=C_n,\end{array}$$ where $C_n=C(n;h,K)$ is a positive and finite constant if ${\int }_{\mathbb{T}}{K}_{x,h}(t)\nu (dx)<\mathrm{\infty }$ for all $t\in \mathbb{T}$, and $\nu$ is a count or Lebesgue measure on $\mathbb{T}$.

It is easy to see that $C_n=1$ for the estimators (6) with DiracDU kernel or any classical (symmetric) associated kernel. Indeed, for the DiracDU kernel estimation we have $$\begin{array}{rl}{\displaystyle \sum _{x=0}^{c-1}{\hat{f}}_n(x)}& ={\displaystyle \sum _{x=0}^{c-1}\{(1-h){\mathbb{1}}_{\left\{x\right\}}(X_1)+{\displaystyle \frac{h}{c-1}}{\mathbb{1}}_{{\mathbb{S}}_c\setminus \left\{x\right\}}(X_1)\}}\\ & =(1-h)+{\displaystyle \frac{h}{c-1}}(c-1)\\ & =1.\end{array}$$ In general we have $C_n\ne 1$ for other discrete and also continuous associated kernels, but it is always close to $1$. In practice, we compute $C_n$ depending on observations before normalizing ${\hat{f}}_n$ to be a p.m.f. or a p.d.f. The following code helps to compute the normalizing constant, e.g., for gamma kernel estimation:

R> data("faithful", package = "datasets")
R> x <- faithful$waiting
R> f <- dke.fun(x, ker = "GA", 0.1)
R> f$C_n
    
[1] 0.9888231

Without loss of generality, we study $x\mapsto {\hat{f}}_n(x)$ up to a normalizing constant which is used at the end of the density estimation process. Notice that that non-classical associated kernel estimators ${\hat{f}}_n$ are improper density estimates or as kind of “balloon estimators”; see (Sain 2002). There are two ways to normalize these estimators (6). The first method is the global normalization using $C_n$ of (7): $$\begin{array}{}\text{(8)}& {\tilde{f}}_n(x)=\frac{{\hat{f}}_n(x)}{{\displaystyle {\int }_{inf(\mathbb{T})}^{sup(\mathbb{T})}{\hat{f}}_n(x)\nu (dx)}},x\in \mathbb{T}.\end{array}$$ Another alternative is to use an adaptive normalization of (6) according to each target $x$: $$\begin{array}{rl}\widetilde{{\stackrel{~}{f}}_n}(x)& ={\displaystyle \frac{1}{n}}{\displaystyle \sum _{i=1}^n\frac{K_{x,h}(X_i)}{{\displaystyle {\int }_{inf(\mathbb{T})}^{sup(\mathbb{T})}{K}_{x,h}(X_i)\nu (dx)}},x\in \mathbb{T},}\end{array}$$ but this approach, with similar results than (8), is not used here. The representations are done with the global normalization (8). In the package, we also compute the normalizing constant (7) for any data set.

In discrete cases, the integrated squared error (ISE) defined by $${ISE}_0=\sum _{x\in \mathbb{N}}\{{\tilde{f}}_n(x)-f_0(x){\}}^2,$$ is the criteria used to measure numerically the discrete smoothness of ${\tilde{f}}_n$ from (8) with the empirical or naive p.m.f. $f_0$ such that $\sum _{x\in \mathbb{N}}{f}_0(x)=1$; see, e.g., (Kokonendji and Senga Kiessé 2011). Concerning the continuous variables, the histogram gives a graphical measure of comparison with ${\tilde{f}}_n$; see, for example, Figure 2.

Some theoretical aspects of the normalizing constant

In this section, we present some theoretical aspects of the normalizing constant $C_n$ of (7) and two examples in the continuous and discrete cases. We first recall the following result on pointwise properties of the estimator (6).

Lemma 2.3.2

(Kokonendji and Senga Kiessé 2011),(Libengué 2013) Let $x\in \mathbb{T}$ be a target and $h\equiv h_n$ a bandwidth. Assuming $f$ is in the class ${\mathcal{C}}^2(\mathbb{T})$ in the continuous case, then $$\begin{array}{}\text{(9)}& Bias\{{\hat{f}}_n(x)\}=A(x,h)f^{\prime }(x)+\frac{1}{2}\{A^2(x,h)+B(x,h)\}{f}^{″}(x)+o(h^2).\end{array}$$ Similar expressions (9) hold in the discrete case, except that $f^{\prime }$ and $f^{″}$ are finite differences of the first and second order respectively.

Furthermore, for the continuous case, if $f$ is bounded on $\mathbb{T}$ then there exists $r_2=r_2\left(K_{x,h}\right)>0$ the largest real number such that ${\Vert K_{x,h}\Vert }_2^2:={\int }_{{\mathbb{S}}_{x,h}}{K}_{x,h}^2(u)du\le c_2(x)h^{-r_2}$, $0\le c_2(x)\le \mathrm{\infty }$ and $$\begin{array}{}\text{(10)}& \mathbb{V}\{{\hat{f}}_n(x)\}=\frac{1}{n}f(x){\Vert K_{x,h}\Vert }_2^2+o\left(\frac{1}{n{h}^{r_2}}\right).\end{array}$$ For discrete situations, the result (10) becomes $$\begin{array}{}& \mathbb{V}\{{\hat{f}}_n(x)\}={\displaystyle \frac{1}{n}}f(x)[\{\mathbb{P}({\mathcal{Z}}_{x,h}=x){\}}^2-f(x)],\end{array}$$ where ${\mathcal{Z}}_{x,h}$ denotes the discrete random variable with p.m.f.  $K_{x,h}$.

It is noticeable that the bias (9) is bigger than the one with symmetric kernels and thus can be reduced; see, e.g., (Zhang 2010), (Zhang and Karunamuni 2010) and (Libengué 2013).

Proposition 2.3.3

Following notations in Lemma 2.3.2, the mean and variance of $C_n$ of (7) are respectively: $$\begin{array}{}\text{(11)}& \mathbb{E}\left(C_n\right)\simeq 1+{\int }_{inf(\mathbb{T})}^{sup(\mathbb{T})}\{A(x,h)f^{\prime }(x)+\frac{1}{2}[A^2(x,h)+B(x,h)]f^{″}(x)\}\nu (dx),\end{array}$$

$$\begin{array}{}\text{(12)}& \mathbb{V}\left(C_n\right)\simeq \{\begin{array}{ll}{\displaystyle \frac{1}{n}}{\displaystyle {\int }_{inf(\mathbb{T})}^{sup(\mathbb{T})}(f(x){\Vert K_{x,h}\Vert }_2^2)dx}& \text{ if }\mathbb{T}\text{ is }\text{ continuous, }\\ \\ {\displaystyle \frac{1}{n}}{\displaystyle \sum _{x\in \mathbb{T}}\left(f(x)[\{\mathbb{P}({\mathcal{Z}}_{x,h}=x){\}}^2-f(x)]\right)}& \text{ if }\mathbb{T}\text{ is }\text{ discrete, }\end{array}\end{array}$$ where $inf(\mathbb{T})$ and $sup(\mathbb{T})$ are respectively the infimum and supremum of $\mathbb{T}$, the measure $\nu$ is Lesbesgue or count on the support $\mathbb{T}$, and where “$\simeq$” stands for approximation.

Proof. From Lemma 2.3.2 and the Fubini theorem, we successively show (11) as follows: $$\begin{array}{rl}\mathbb{E}\left(C_n\right)& =\mathbb{E}\left({\displaystyle {\int }_{inf(\mathbb{T})}^{sup(\mathbb{T})}{\hat{f}}_n(x)\nu (dx)}\right)={\displaystyle {\int }_{inf(\mathbb{T})}^{sup(\mathbb{T})}\mathbb{E}\left({\hat{f}}_n(x)\right)\nu (dx)}\\ & ={\displaystyle {\int }_{inf(\mathbb{T})}^{sup(\mathbb{T})}(Bias\{{\hat{f}}_n(x)\}+f(x))\nu (dx)}\\ & \simeq {\int }_{inf(\mathbb{T})}^{sup(\mathbb{T})}\{A(x,h)f^{\prime }(x)+\frac{1}{2}[A^2(x,h)+B(x,h)]f^{″}(x)+f(x)\}\nu (dx)\\ & \simeq {\displaystyle {\int }_{inf(\mathbb{T})}^{sup(\mathbb{T})}f(x)\nu (dx)+{\displaystyle {\int }_{inf(\mathbb{T})}^{sup(\mathbb{T})}\{A(x,h)f^{\prime }(x)+\frac{1}{2}[A^2(x,h)+B(x,h)]f^{″}(x)\}\nu (dx)}}\\ & \simeq 1+{\int }_{inf(\mathbb{T})}^{sup(\mathbb{T})}\{A(x,h)f^{\prime }(x)+\frac{1}{2}[A^2(x,h)+B(x,h)]f^{″}(x)\}\nu (dx).\end{array}$$

The variance (12) is trivial from Lemma 2.3.2. $◻$

Example 1. Let $f$ be an exponential density with parameter $\gamma >0$. Thus, one has: $$\begin{array}{rl}f(x)& =\gamma \exp (-\gamma x),\\ f^{\prime }(x)& =-{\gamma }^2\exp (-\gamma x),\\ f^{″}(x)& ={\gamma }^3\exp (-\gamma x).\end{array}$$ Consider the lognormal kernel with $A(x,h)=x(\exp (3h^2/2)-1)$, $B(x,h)=x^2\exp (3h^2)(\exp (h^2)-1)$ and ${\Vert L{N}_{x,h}\Vert }_2^2=1/(2\pi h{x}^{1/2})$; see (Libengué 2013). Then, using the Taylor formula around $h$, the expressions of $\mathbb{E}\left(C_n\right)$ and $\mathbb{V}\left(C_n\right)$ are: $$\mathbb{E}\left(C_n\right)\simeq 1+[(-1+\exp \frac{3h^2}{2})+\frac{1}{2}\{{(-1+\exp \frac{3h^2}{2})}^2+\exp 3h^2(-1+\exp h^2)\}]\simeq 1-\frac{h^2}{2}$$ and $\mathbb{V}\left(C_n\right)\simeq \gamma (2nh\sqrt{\pi }{)}^{-1}{\int }_0^{\mathrm{\infty }}{z}^{-1}\exp (-z)dz$ with ${\int }_0^{\mathrm{\infty }}{z}^{-1}\exp (-z)dz\approx 16.2340$ by computation with R. Thus, the quantity $C_n$ cannot be equal to 1.

Example 2. Let $f$ be a Poisson p.m.f. with parameter $\lambda$ and thus $f(x)={\lambda }^x\exp (-\lambda )/x!$. The finite differences $f^{(k)}(x)$ of order $k\in \{1,2,\dots \}$ at $x\in \mathbb{N}$ are given by the recursive relation: $$f^{(k)}(x)=\{f^{(k-1)}(x){\}}^{(1)}\text{ with }{f}^{(1)}(x)=\{\begin{array}{ll}\{f(x+1)-f(x-1)\}/2,\text{ if }x\in \mathbb{N}\setminus \left\{0\right\},& \\ f(1)-f(0),\text{ if }x=0,& \end{array}$$ and $$f^{(2)}(x)=\{\begin{array}{ll}\{f(x+2)-2f(x)+f(x-2)\}/4,\text{ if }x\in \mathbb{N}\setminus \{0,1\},& \\ \{f(3)-3f(1)+f(0)\}/4,\text{ if }x=1,& \\ \{f(2)-2f(1)+f(0)\}/2,\text{ if }x=0.& \end{array}$$

Considering the binomial kernel with $A(x,h)=x+h$, $B(x,h)=(x+h)(1-h)/(x+1)$, we successively obtain $$\begin{array}{rl}\mathbb{E}\left(C_n\right)& \simeq 1+\frac{h}{4}(f(3)+3f(2)-3f(1)-2f(0)\\ & \phantom{\simeq }+{\displaystyle \sum _{x=2}^{\mathrm{\infty }}2\{f(x+1)-f(x-1)\}+{\displaystyle \frac{(x+3)\{f(x+2)-2f(x)+f(x-2)\}}{x+1}})}\\ & \phantom{\simeq }+(\frac{3f(3)+8f(2)-17f(1)+3f(0)}{2}\\ & \phantom{\simeq }+{\displaystyle \sum _{x=2}^{\mathrm{\infty }}\frac{x\{f(x+1)-f(x-1)\}}{2}+{\displaystyle \frac{(x^3+x^2+1)\{f(x+2)-2f(x)+f(x-2)\}}{x+1}}})\\ & \simeq 1+\frac{h\exp (-\lambda )}{4}[\frac{{\lambda }^3}{3!}+3\frac{{\lambda }^2}{2!}-3\lambda -2\\ & \phantom{\simeq }+{\displaystyle \sum _{x=2}^{\mathrm{\infty }}\{\frac{2{\lambda }^{x+1}}{(x+1)!}-\frac{2{\lambda }^{x-1}}{(x-1)!}+{\displaystyle \frac{x+3}{x+1}}(\frac{{\lambda }^{x+2}}{(x+2)!}-2\frac{{\lambda }^x}{x!}+\frac{{\lambda }^{x-2}}{(x-2)!})\}}]\\ & \phantom{\simeq }+[\frac{{\lambda }^3}{4}+2{\lambda }^2-\frac{17\lambda }{2}+3\\ & \phantom{\simeq }+{\displaystyle \sum _{x=2}^{\mathrm{\infty }}\{\frac{x}{2}(\frac{{\lambda }^{x+1}}{(x+1)!}-\frac{{\lambda }^{x-1}}{(x-1)!})+{\displaystyle \frac{x^3+x^2+1}{x+1}}(\frac{{\lambda }^{x+2}}{(x+2)!}-2\frac{{\lambda }^x}{x!}+\frac{{\lambda }^{x-2}}{(x-2)!})\}}]\end{array}$$ and $$\mathbb{V}\left(C_n\right)\simeq {\displaystyle \frac{\exp (-\lambda )}{n}}{\displaystyle \sum _{x\in \mathbb{T}}\left(\frac{{\lambda }^x}{x!}[{\{(1+h){\left(\frac{x+h}{x+1}\right)}^x\}}^2-\frac{{\lambda }^x\exp (-\lambda )}{x!}]\right).}$$

Bandwidth selection

Now, we consider the bandwidth selection problems which are generally crucial in nonparametric estimation. Several methods already existing for continuous kernels can be adapted to the discrete case as the classical least-squares cross-validation method; see, for example, (Bowman 1984), (Marron 1987) and references therein. Here, we simply propose three procedures for the bandwidth selection: cross-validation, Bayesian local for binomial and adaptive for the gamma kernel. Also, a review of bayesian bandwidth selection methods is presented. Each time, the smoothing parameter selection is done with the non-normalized version ${\hat{f}}_n$ of the estimator (6) before the global normalization ${\tilde{f}}_n$ of (8).

Cross-validation for any associated kernel

For a given associated kernel $K_{x,h}$ with $x\in \mathbb{T}$ and $h>0$, the optimal bandwidth $h_{cv}$ of $h$ is obtained by cross-validation as $h_{cv}=\arg \underset{\mathit{h}>0}{min}CV(h)$ with $$CV(h)={\displaystyle {\int }_{x\in \mathbb{T}}{\{{\hat{f}}_n(x)\}}^2\nu (dx)-\frac{2}{n}{\displaystyle \sum _{i=1}^n{\hat{f}}_{n,-i}(X_i),}}$$ where ${\hat{f}}_{n,-i}(X_i)=(n-1{)}^{-1}{\displaystyle \sum _{j\ne i}{K}_{X_i,h}(X_j)}$ is being computed as ${\hat{f}}_n(X_i)$ by excluding the observation $X_i$ and $\nu$ is the Lebesgue or count measure. This method is applied to all estimators (6) with associated kernels cited in this paper, independently on the support $\mathbb{T}$ of $f$ to be estimated.

Table 2: Summary of arguments and results of hcvc.fun.

Arguments	Description
Vec	The positive continuous data sample.
seq.bws	The sequence of bandwidths where to compute the cross-validation function.
ker	The associated kernel.
a0,a1	The bounds of the support of extended beta kernel. Default values are
	respectively $0$ and $1$.
a	The arm of the discrete triangular kernel. Default value is $1$.
c	The number of categories in DiracDU kernel. Default value is $2$.
Results	Description
hcv	The optimal bandwidth obtained by cross-validation.
seq.h	The sequence of bandwidths used to compute hcv.
CV	The values of the cross-validation function.

Table 2 gives the arguments and results of the cross-validation function hcvc.fun defined for continuous data are below. The hcvd.fun is the corresponding function for discrete data. The hcvc.fun is performed with the Old Faithful geyser data described in (Azzalini and Bowman 1990) and (Härdle 2012). The dataset concerns waiting time between eruptions and the duration of the eruption for the Old Faithful geyser in Yellowstone National Park, Wyoming, USA. The following codes and Figure 2 give smoothing density estimation with various associated kernels of the waiting time variable.

R> data("faithful", package = "datasets")
R> x <- faithful$waiting
R> f1 <- dke.fun(x, 0.1, "continuous", ker = "GA")
R> f2 <- dke.fun(x, 0.036, "continuous", ker = "LN") 
R> f3 <- dke.fun(x, 0.098, "continuous", ker = "RIG")
R> f4 <- dke.fun(x, 0.01, "continuous", ker = "BE", a0 = 40, a1 = 100)
R> t <- seq(min(x), max(x), length.out = 100)
R> hist(x, probability = TRUE, xlab = "Waiting times (in min.)",
+    ylab = "Frequency", main = "", border = "gray")
R> lines(t, f1$fn, lty = 2, lwd = 2, col = "blue")
R> lines(t, f2$fn, lty = 5, lwd = 2, col = "black")
R> lines(t, f3$fn, lty = 1, lwd = 2, col = "green")
R> lines(t, f4$fn, lty = 4, lwd = 2, col = "grey")
R> lines(density(x, width = 12), lty = 8, lwd = 2, col = "red")
R> legend("topleft", c("Gamma", "Lognormal", "Reciprocal inverse Gaussian", 
+    "Extended beta", "Gaussian"), col = c("blue", "black", "green", "grey", "red"),
+    lwd = 2, lty = c(2, 5, 1, 4, 8), inset = .0)

Figure 2: Smoothing density estimation of the Old Faithful geyser data (Azzalini and Bowman 1990) by some continuous associated kernels with the support of observations $[43,96]=\mathbb{T}$.

A review of Bayesian bandwidth selection

Bayesian inference grows out of the simple formula known as Bayes rule. Assume we have two random variables $A$ and $B$. A principle rule of probability theory known as the chain rule allows us to specify the joint probability of $A$ and $B$ taking on particular values $a$ and $b$, $P(a,b)$, as the product of the conditional probability that $A$ will take on value $a$ given that $B$ has taken on value $b$, $P(a|b)$, and the marginal probability that $B$ takes on value $b$, $P(b)$. Which gives us:
Joint probability = Conditional Probability $\times$ Marginal Probability.
Thus we have: $P(a,b)=P(a|b)P(b).$
This expression (Bayes rule) indicates that we can compute the conditional probability of a variable $A$ given the variable $B$ from the conditional probability of $B$ given $A$. This introduces the notion of prior and posterior knowledge.

Prior and posterior knowledge.

A prior probability is the probability available to us beforehand, and before making any additional observations. A posterior probability is the probability obtained from the prior probability after making additional observation to the prior knowledge available.

Summarizing the Bayesian approach.

The Bayesian approach to parameter estimation works as follows:

1. Formulate our knowledge about a situation.

2. Gather data.

3. Obtain posterior knowledge that updates our beliefs.

How do we formulate our knowledge about a situation?

1. Define a distribution model which expresses qualitative aspects of our knowledge about the situation. This model will have some unknown parameters, which will be dealt with as random variables.

2. Specify a prior probability distribution which expresses our subjective beliefs and subjective uncertainty about the unknown parameters, before seeing the data.

After gathering the data, how do we obtain posterior knowledge?

1. Compute posterior probability distribution which estimates the unknown parameters using the rules of probability and given the observed data, presenting us with updated beliefs.

(Zougab et al. 2013) proposed a Bayesian approach based upon a likelihood cross-validation approximation and a Markov chain Monte Carlo (MCMC) method for deriving the global optimal bandwidth using the famous binomial kernel. However, a global bandwidth does not generally provide a good estimator for complex p.m.f.’s, in particular for small and moderate sample sizes. Generally, the global discrete associated kernel estimator tends to simultaneously under- and oversmooth $f(x)$.

In order to improve the global discrete associated kernel estimator, in particular for complex count data with small and moderate sample sizes, (Zougab et al. 2012) and (Zougab et al. 2013) adapted two versions of variable bandwidths for discrete associated kernel estimator and proposed Bayesian approaches for selecting these variable bandwidths. Note that these two versions are originally proposed for kernel density estimation (see, e.g., Sain and Scott (1996); Breiman et al. (1977); Abramson (1982); Brewer (2000); Zhang et al. (2006); Zougab et al. (2014b) and Zhang et al. (2016)).

Recently, (Zougab et al. 2012) have considered the local discrete associated kernel estimator (balloon estimator in discrete case) and have derived the closed form of the variable bandwidth at each point $x$ for which the p.m.f. is estimated by considering the binomial kernel estimator and locally treating the bandwidth as a random quantity with a beta prior distribution. This approach outperforms existing classical global methods, namely, MISE and CV in particular for small and moderate sample sizes. (Zougab et al. 2013) have also proposed the adaptive discrete associated kernel estimator (sample-point estimator in discrete situation), which replaces $h$ by $h_i$ for each observation $x_i$ with $i=1,\dots ,n$, and then employs the Bayesian approach for estimating the adaptive bandwidths $h_i$. The authors have considered the binomial kernel and the beta prior for each variable bandwidth $h_i$, and have shown that this approach performs better than the popular classical global selectors.

Bayesian estimation of localized bandwidth for the binomial kernel

An alternative to the cross-validation for bandwidth selection is by using Bayesian methods. These methods have been investigated with three different procedures: local, global and adaptive; see respectively (Zougab et al. 2012),(Zougab et al. 2013),(Zougab et al. 2014a). In terms of integrated squared error and execution times, the local Bayesian outperforms the other Bayesian procedures. In the local Bayesian framework, the variable bandwidth is treated as parameter with prior $\pi (\cdot )$. Under squared error loss function, the Bayesian bandwidth selector is the posterior mean; see (Zougab et al. 2012).

First, as we have mentioned above, $f(x)$ can be approximated by $$f(x|h)=f_h(x)={\displaystyle \sum _{u\in \mathbb{T}}f(u)B_{x,h}(u)=\mathbb{E}\{B_{x,h}(X)\},}$$ where $B_{x,h}$ is the binomial kernel and $X$ is a random variable with p.m.f. $f$. Now, considering $h$ as a scale parameter for $f_h(x)$, the local approach consists of using $f_h(x)$ and constructing a Bayesian estimator for $h(x)$.

Indeed, let $\pi (h)$ denote the beta prior density of $h$ with positive parameters $\alpha$ and $\beta$. By the Bayes theorem, the posterior of $h$ at the point of estimation $x$ takes the form $$\pi (h|x)={\displaystyle \frac{f_h(x)\pi (h)}{\int f_h(x)\pi (h)dh}}.$$ Since $f_h$ is unknown, we use ${\hat{f}}_h$ as natural estimator of $f_h$, and hence we can estimate the posterior by $$\pi (h|x,X_1,X_2,\dots ,X_n)={\displaystyle \frac{{\hat{f}}_h(x)\pi (h)}{\int {\hat{f}}_h(x)\pi (h)dh}}.$$ Under the squared error loss, the Bayes estimator of the smoothing parameter $h(x)$ is the posterior mean and is given by ${\hat{h}}_n(x)=\int h\hat{\pi }(h|x,X_1,X_2,\dots ,X_n)dh$. Exact approximation is $${\hat{h}}_n(x)={\displaystyle \frac{{\displaystyle \sum _{i=0}^n\sum _{k=0}^{X_i}\frac{x^k}{(x+1-X_i)!k!(X_i-k)!}B(X_i+\alpha -k+1,x+\beta +1-X_i)}}{{\displaystyle \sum _{i=0}^n\sum _{k=0}^{X_i}\frac{x^k}{(x+1-X_i)!k!(X_i-k)!}B(X_i+\alpha -k,x+\beta +1-X_i)}}},\mathrm{\forall }x\in \mathbb{N}\text{ with }{X}_i\le x+1,$$ where $B(\cdot ,\cdot )$ is the beta function; see (Zougab et al. 2012) for more details.

Bayesian estimation of adaptive bandwidth for the gamma kernel

The bandwidth $h$ in the gamma kernel density estimation can be allowed to be adaptive. This approach gives a variable bandwidth $h_i$ for each observation $X_i$ in place of the initial fixed bandwidth $h$. Following (Zougab et al. 2014b), we suggest using Bayesian methods to estimate such adaptive bandwidths or variable bandwidths $h_i,i=1,\dots ,n$. Thus, we treat $h_i$ as a random variable with a prior distribution $\pi (\cdot )$. The estimator ($\text{6}$) with gamma kernel of Section (5) and variable bandwidths are reformulated as $$\begin{array}{}\text{(13)}& {\hat{f}}_n(x)=\frac{1}{n}\sum _{i=1}^{n}G{A}_{x,h_i}(X_i).\end{array}$$ The leave-one-out kernel estimator of $f(X_i)$ deduced from ($\text{13}$) is $$\begin{array}{}\text{(14)}& \hat{f}(X_i\mid \{X_{-i}\},h_i)=\frac{1}{n-1}\sum _{j=1,j\ne i}^{n}G{A}_{X_i,h_i}(X_j),\end{array}$$ where $\{X_{-i}\}$ denotes the set of observations excluding $X_i$. The posterior distribution for each variable bandwidth $h_i$ given $X_i$ provided from the Bayesian rule is expressed as follow $$\begin{array}{}\text{(15)}& \pi (h_i\mid X_i)=\frac{\hat{f}(X_i\mid \{X_{-i}\},h_i)\pi (h_i)}{{\int }_0^{\mathrm{\infty }}\hat{f}(X_i\mid \{X_{-i}\},h_i)\pi (h_i)d{h}_i}.\end{array}$$ We obtain the Bayesian estimator ${\tilde{h}}_i$ of $h_i$ by using the quadratic loss function $$\begin{array}{}\text{(16)}& {\tilde{h}}_i=\mathbb{E}(h_i\mid X_i).\end{array}$$ In the following, we assume that each $h_i=h_i(n)$ has an inverse gamma prior distribution $\mathcal{IGA}(\alpha ,\beta )$ with the shape parameter $\alpha >0$ and scale parameter ${\beta }_{}>0$. The density of $\mathcal{IGA}(a,b)$ with $a,b>0$ is defined as $$\begin{array}{}\text{(17)}& {\mathrm{\Phi }}_{a,b}(z)=\frac{b^a}{\mathrm{\Gamma }(a)}{z}^{-a-1}\exp (-b/z){\mathbb{1}}_{(0,\mathrm{\infty })}(z).\end{array}$$ This allows us to obtain the closed form of the posterior density and the Bayesian estimator given by the following result.

Theorem 2.3.4

For fixed $i\in \{1,2,\dots ,n\}$, consider each observation $X_i$ with its corresponding bandwidth $h_i$. Using the gamma kernel estimator ((13)) and the inverse gamma prior distribution $\mathcal{IGA}(\alpha ,\beta )$ given in ((17)) with $\alpha >1/2$ and $\beta >0$ for each $h_i$, then:

1. The posterior density ($\text{15}$) is the following weighted sum of inverse gamma $$\pi (h_i\mid X_i)=\frac{1}{D_{ij}}\sum _{j=1,j\ne i}^n\{A_{ij}{\mathrm{\Phi }}_{\alpha +1/2,B_{ij}}(h_i){\mathbb{1}}_{(0,\mathrm{\infty })}(X_i)+C_j{\mathrm{\Phi }}_{\alpha +1,X_j+\beta }(h_i){\mathbb{1}}_{\{0\}}(X_i)\}$$ with $A_{ij}=[\mathrm{\Gamma }(\alpha +1/2)]/({\beta }^{\alpha }{X}_i^{1/2}\sqrt{2\pi }{B}_{ij}^{\alpha +1/2})$, $B_{ij}=X_i\log X_i-X_i\log X_j+X_j-X_i+\beta$, $C_j=\mathrm{\Gamma }(\alpha +1)/[{\beta }_{}^{-\alpha }(X_j+\beta {)}^{\alpha +1}]$ and $D_{ij}=\sum _{j=1,j\ne i}^n\{A_{ij}{\mathbb{1}}_{(0,\mathrm{\infty })}(X_i)+C_j{\mathbb{1}}_{\{0\}}(X_i)\}$.

2. The Bayesian estimator ${\tilde{h}}_i$ of $h_i$, given in ($\text{16}$), is $${\tilde{h}}_i=\frac{1}{D_{ij}}\sum _{j=1,j\ne i}^n\{\frac{A_{ij}{B}_{ij}}{\alpha -1/2}{\mathbb{1}}_{(0,\mathrm{\infty })}(X_i)+\frac{(X_j+\beta )C_j}{\alpha }{\mathbb{1}}_{\{0\}}(X_i)\}$$ according to the previous notations of $A_{ij}$, $B_{ij}$, $C_j$ and $D_{ij}$.

Proof. (i) Let us represent $\pi (h_i\mid X_i)$ of ((15)) as the ratio of $N(h_i\mid X_i):=\hat{f}(X_i\mid \{X_{-i}\},h_i)\pi (h_i)$ and ${\int }_0^{\mathrm{\infty }}N(h_i\mid X_i)d{h}_i$. From ((14)) and ((17)) the numerator is, first, equal to $$\begin{array}{}\text{(18)}& \begin{array}{rl}N(h_i\mid X_i)& =(\frac{1}{n-1}\sum _{j=1,j\ne i}^{n}G{A}_{X_i,h_i}(X_j))(\frac{{\beta }_{}^{\alpha }}{\mathrm{\Gamma }(\alpha )}{h}_i^{-\alpha -1}\exp (-\beta /h_i))\\ & =\frac{[\mathrm{\Gamma }(\alpha ){]}^{-1}}{(n-1)}\sum _{j=1,j\ne i}^n\frac{G{A}_{X_i,h_i}(X_j)}{{\beta }_{}^{-\alpha }{h}_i^{\alpha +1}}\exp (-{\beta }_{}/h_i).\end{array}\end{array}$$

Following (Chen 2000), we assume that for all $X_i\in (0,\mathrm{\infty })$ one has $1+(X_i/h_i)\to \mathrm{\infty }\text{ as }n\to \mathrm{\infty }$. Using the Stirling formula $\mathrm{\Gamma }(z+1)\simeq \sqrt{2\pi }{z}^{z+1/2}\exp (-z)\text{ as }z\to \mathrm{\infty }$, the term of the sum in ((18)) can be successively calculated as

$$\begin{array}{}\text{(19)}& \begin{array}{rl}\frac{G{A}_{X_i,h_i}(X_j)}{{\beta }_{}^{-\alpha }{h}_i^{\alpha +1}}\exp (-{\beta }_{}/h_i)& =\frac{X_j^{(X_i/h_i)}\exp (-X_j/h_i)}{h_i^{1+(X_i/h_i)}\mathrm{\Gamma }[1+(X_i/h_i)]{\beta }^{-\alpha }{h}_i^{\alpha +1}}\exp (-\beta /h_i)\\ & =\frac{\exp [-(X_j+{\beta }_{}-X_i\log X_j)/h_i]}{{\beta }^{-\alpha }{h}_i^{(X_i/h_i)+\alpha +2}\sqrt{2\pi }\exp (-X_i/h_i)(X_i/h_i{)}^{(X_i/h_i)+1/2}}\\ & =\frac{\mathrm{\Gamma }(\alpha +1/2)}{{\beta }_{}^{-\alpha }{X}_i^{1/2}\sqrt{2\pi }{B}_{ij}^{\alpha +1/2}}\times \frac{B_{ij}^{\alpha +1/2}\exp [-B_{ij}/h_i]}{h_i^{\alpha +3/2}\mathrm{\Gamma }(\alpha +1/2)}\\ & =A_{ij}{\mathrm{\Phi }}_{\alpha +1/2,B_{ij}}(h_i),\end{array}\end{array}$$

with $B_{ij}=X_i\log X_i-X_i\log X_j+X_j-X_i+{\beta }_{}$, $A_{ij}=[X_j^{-1}\mathrm{\Gamma }(\alpha +1/2)]/({\beta }_{}^{-\alpha }{X}_i^{-1/2}\sqrt{2\pi }{B}_{ij}^{\alpha +1/2})$ and ${\mathrm{\Phi }}_{\alpha +1/2,B_{ij}}(h_i)$ is given in ((17)).

Also, for $X_i=0$, the term of the sum ((18)) can be expressed as follows

$$\begin{array}{}\text{(20)}& \begin{array}{rl}\frac{G{A}_{0,h_i}(X_j)}{{\beta }_{}^{-\alpha }{h}_i^{\alpha +1}}\exp (-{\beta }_{}/h_i)& =\frac{\exp (-X_j/h_i)}{{\beta }^{-\alpha }{h}_i^{\alpha +2}}\exp (-\beta /h_i)\\ & =\frac{\mathrm{\Gamma }(\alpha +1)}{{\beta }_{}^{-\alpha }(X_j+\beta {)}^{\alpha +1}}\times \frac{(X_j+\beta {)}^{\alpha +1}\exp [-(X_j+\beta )/h_i]}{h_i^{\alpha +2}\mathrm{\Gamma }(\alpha +1)}\\ & =C_j{\mathrm{\Phi }}_{\alpha +1,X_j+\beta }(h_i),\end{array}\end{array}$$

with $C_j=\mathrm{\Gamma }(\alpha +1)/[{\beta }_{}^{-\alpha }(X_j+\beta {)}^{\alpha +1}]$ and ${\mathrm{\Phi }}_{\alpha +1,X_j+\beta }(h_i)$ is given in ((17)). Combining (19) and (20), the expression of $N(h_i\mid X_i)$ in ((18)) becomes

$$\begin{array}{}\text{(21)}& \begin{array}{r}N(h_i\mid X_i)=\frac{[\mathrm{\Gamma }(\alpha ){]}^{-1}}{(n-1)}\sum _{j=1,j\ne i}^n\{A_{ij}{\mathrm{\Phi }}_{\alpha +1/2,B_{ij}}(h_i){\mathbb{1}}_{(0,\mathrm{\infty })}(X_i)+C_j{\mathrm{\Phi }}_{\alpha +1,X_j+\beta }(h_i){\mathbb{1}}_{\{0\}}(X_i)\}.\end{array}\end{array}$$

From (21), the denominator is successively computed as follows:

$$\begin{array}{}\text{(22)}& \begin{array}{rl}{\int }_0^{\mathrm{\infty }}N(h_i\mid X_i)d{h}_i& =\frac{[\mathrm{\Gamma }(\alpha ){]}^{-1}}{(n-1)}\sum _{j=1,j\ne i}^n(A_{ij}{\int }_0^{\mathrm{\infty }}{\mathrm{\Phi }}_{\alpha +1/2,B_{ij}}(h_i){\mathbb{1}}_{(0,\mathrm{\infty })}(X_i)d{h}_i\\ & \phantom{=}+C_j{\int }_0^{\mathrm{\infty }}{\mathrm{\Phi }}_{\alpha +1,X_j+\beta }(h_i){\mathbb{1}}_{\{0\}}(X_i)d{h}_i)\\ & =\frac{[\mathrm{\Gamma }(\alpha ){]}^{-1}}{(n-1)}\sum _{j=1,j\ne i}^n\{A_{ij}{\mathbb{1}}_{(0,\mathrm{\infty })}(X_i)+C_j{\mathbb{1}}_{\{0\}}(X_i)\}\\ & =\frac{[\mathrm{\Gamma }(\alpha ){]}^{-1}}{(n-1)}{D}_{ij},\end{array}\end{array}$$

with $D_{ij}=\sum _{j=1,j\ne i}^n(A_{ij}{\mathbb{1}}_{(0,\mathrm{\infty })}(X_i)+C_j{\mathbb{1}}_{\{0\}}(X_i))$. Finally, the ratio of ((21)) and ((22)) leads to the result of Part (i).

(ii) Let us remember that the mean of the inverse gamma distribution $\mathcal{IG}(\alpha ,\beta )$ is $\beta /(\alpha -1)$. Thus, the expression of $\pi (h_i\mid X_i)$ in (15) is given by

$$\pi (h_i\mid X_i)=\frac{1}{D_{ij}}\sum _{j=1,j\ne i}^n\{A_{ij}{\mathrm{\Phi }}_{\alpha +1/2,B_{ij}}(h_i){\mathbb{1}}_{(0,\mathrm{\infty })}(X_i)+C_j{\mathrm{\Phi }}_{\alpha +1,X_j+\beta }(h_i){\mathbb{1}}_{\{0\}}(X_i)\}$$ and, therefore, ${\tilde{h}}_i=\mathbb{E}(h_i\mid X_i)={\int }_0^{\mathrm{\infty }}{h}_i\pi (h_i\mid X_i)d{h}_i$ is finally $${\tilde{h}}_i=\mathbb{E}(h_i\mid X_i)=\frac{1}{D_{ij}}\sum _{j=1,j\ne i}^n\{\frac{A_{ij}{B}_{ij}}{\alpha -1/2}{\mathbb{1}}_{(0,\mathrm{\infty })}(X_i)+\frac{(X_j+\beta )C_j}{\alpha }{\mathbb{1}}_{\{0\}}(X_i)\}.$$

This new method of selecting bandwidth by the Bayesian adaptive procedure will be implemented in a future version of the Ake package.

4 Bandwidth selection for kernel regression involving associated kernels

One of the most often encountered models in nonparametric statistics is the regression model. The function that provides the best prediction of a dependent variable $y$ in terms of an independent variable $x$ is the conditional expectation $\mathbb{E}(y/x)=m(x)$. This is called regression function and its estimation from a sequence of $n$ pairs $(x_i,y_i)$, $i=1,\dots ,n$ is a problem in statistics. We will consider the case $(x,y)\in \mathbb{T}\times \mathbb{R}$. For simplicity, we take $\mathbb{T}=\mathbb{R}$ if $x$ is a continuous variable and $\mathbb{T}=\mathbb{N}$ in the discrete case. The classical non parametric regression model between two variables $y$ and $x$ is $$\begin{array}{}\text{(23)}& y_i=m(x_i)+{ϵ}_i,\end{array}$$ where $\mathbf{\text{y}}=(y_1,\dots ,y_n)$ is a response vector, $\mathbf{\text{x}}=(x_1,\dots ,x_n)$ is an explanatory vector, $ϵ=({ϵ}_1,\dots ,{ϵ}_n)$ is the error following a Gaussian distribution with zero mean and finite variance ${\sigma }^2$, i.e., ${ϵ}_i\sim \mathcal{N}(0,{\sigma }^2)$ and $m:\mathbb{T}\mapsto \mathbb{R}$ is the unknown regression function. Several methods have been proposed to estimate the regression function in the continuous case. We cite for example the histograms introduced by (Tukey 1961) and studied by (Geoffroy 1980) and (Lecoutre 1990), the spline method which can be found in (Reinsch 1967), (Silverman 1985) and (Wahba 1990) and also the regression using partition proposed by (Breiman et al. 1984).

As for the density or probability mass function, the estimate of the regression function by the kernel method is the most used because of its good asymptotic properties and interest in practice. Introduced initially for continuous density estimation by (Rosenblatt 1956) and (Parzen 1962), this method was adopted by (Nadaraya 1964) and (Watson 1964) for estimating the continuous regression function. It was also applied to smooth the discrete regression function $m$ for $x\in \mathbb{N}$. Some studies have been done to estimate the discrete regression function, using the Dirac-type kernel (naive estimator) or discrete kernels of (Aitchison and Aitken 1976). However, the naive estimator is appropriate only when the sample size is large, and the discrete kernel of (Aitchison and Aitken 1976) is only suitable for categorical data; see (Hayfield and Racine 2008) and also (Hayfield and Racine 2014). (Kokonendji et al. 2009) adapted the (Nadaraya 1964) and (Watson 1964) kernel to the discrete unknown function $m$, using the discrete associated kernels. In their work, using the integrated mean square error and the coefficient of determination $R^2$, they showed that the binomial or discrete triangular kernels are better compared to the optimal Epanechnikov kernel. In this section we present the theoretical foundations of the estimated regression function with continuous and discrete associated kernels.

Both in continuous and discrete cases, consider the relation between a response variable $Y$ and an explanatory variable $x$ given by $$\begin{array}{}\text{(24)}& Y=m\left(x\right)+ϵ,\end{array}$$ where $m$ is an unknown regression function from $\mathbb{T}\subseteq \mathbb{R}$ to $\mathbb{R}$ and $ϵ$ the disturbance term with null mean and finite variance. Let $(X_1,Y_1),\dots ,(X_n,Y_n)$ be a sequence of i.i.d. random variables on $\mathbb{T}\times \mathbb{R}(\subseteq {\mathbb{R}}^2)$ with $m(x)=\mathbb{E}(Y|X=x)$ of ((24)). Using (continuous or discrete) associated kernels, the (Nadaraya 1964) and (Watson 1964) estimator ${\hat{m}}_n$ of $m$ is $$\begin{array}{}\text{(25)}& {\hat{m}}_n(x;h)={\displaystyle \sum _{i=1}^n{\displaystyle \frac{Y_i{K}_{x,h}\left(X_i\right)}{\sum _{i=1}^n{K}_{x,h}\left(X_i\right)}}={\hat{m}}_n(x),\mathrm{\forall }x\in \mathbb{T}\subseteq \mathbb{R},}\end{array}$$ where $h\equiv h_n$ is the smoothing parameter such that $h_n\to 0$ as $n\to \mathrm{\infty }$.

Besides the criterion of kernel support, we retain the root mean squared error (RMSE) and also the practical coefficient of determination given respectively by $$RMSE=\sqrt{\frac{1}{n}{\displaystyle \sum _{i=1}^n{\{Y_i-{\hat{m}}_n(X_i)\}}^2}}$$ and $$R^2=\frac{\sum _{i=1}^n{\{{\hat{m}}_n(X_i)-\bar{y}\}}^2}{\sum _{i=1}^n(Y_i-\bar{y}{)}^2},$$ with $\bar{y}=n^{-1}(Y_1+\dots +Y_n)$.

In discrete cases, the reg.fun function for (25) is used with the binomial kernel on milk data as follows. This dataset is about average daily fat (kg/day) yields from milk of a single cow for each of the first 35 weeks.

R> data("milk", package = "Ake")
R> x <- milk$week
R> y <- milk$yield
R> h <- reg.fun(x, y, "discrete", "bino", 0.1)
R> h

Bandwidth h:0.1        Coef_det=0.9726

Number of points: 35;   Kernel =  Binomial

data            y         
Min.   : 1.0   Min.   :0.0100  
1st Qu.: 9.5   1st Qu.:0.2750  
Median :18.0   Median :0.3600  
Mean   :18.0   Mean   :0.3986  
3rd Qu.:26.5   3rd Qu.:0.6150  
Max.   :35.0   Max.   :0.7200  
eval.points        m_n         
Min.   : 1.0   Min.   :0.01542  
1st Qu.: 9.5   1st Qu.:0.27681  
Median :18.0   Median :0.35065  
Mean   :18.0   Mean   :0.39777  
3rd Qu.:26.5   3rd Qu.:0.60942  
Max.   :35.0   Max.   :0.70064  

The above reg.fun is also used for continuous cases; see Figure 3 and Table 4 for the motorcycle impact data of (Silverman 1985).

Bandwidth selection

We present two bandwidth selection methods for the regression: the well-known cross-validation for any associated kernel and the Bayesian global for the binomial kernel.

Cross-validation for any associated kernel

For a given associated kernel, the optimal bandwidth parameter is ${\hat{h}}_{cv}={\displaystyle \text{arg }\underset{h>0}{min}\text{ LSCV}(h)}$ with $$\begin{array}{}\text{(26)}& \text{LSCV}(h)=\frac{1}{n}{\displaystyle \sum _{i=1}^n{\{Y_i-{\hat{m}}_{-i}(X_i)\}}^2,}\end{array}$$ where ${\tilde{m}}_{-i}(X_i)$ is computed as ${\hat{m}}_n$ of ((25)) excluding $X_i$; see, e.g., (Kokonendji et al. 2009). The hcvreg.fun function to compute this optimal bandwidth is described in Table 3.

Table 3: Summary of arguments and results of hcvreg.fun.

Arguments	Description
Vec	The explanatory data sample can be discrete or continuous.
y	The response variable.
ker	The associated kernel.
h	The sequence of bandwidths where to compute the optimal bandwidth.
a0,a1	The bounds of the support of extended beta kernel. Default values are
	respectively $0$ and $1$.
a	The arm of the discrete triangular kernel. Default value is $1$.
c	The number of categories in DiracDU kernel. Default value is $2$.
Results	Description
kernel	The associated kernel.
hcv	The optimal bandwidth obtained by cross-validation.
CV	The values of the cross-validation function.
seqbws	The sequence of bandwidths used to compute hcv.

The following code helps to compute the bandwidth parameter by cross-validation on milk data. The associated kernel used is the discrete triangular kernel with arm $a=1$.

R> data("milk", package = "Ake")
R> x <- milk$week
R> y <- milk$yield
R> f <- hcvreg.fun(x, y, type_data = "discrete", ker = "triang", a = 1)
R> f$hcv

[1] 1.141073

When we consider the continuous associated kernel, one needs to set the type of data parameter to “continuous” in the hcvreg.fun function. Thus, the hcvreg.fun and reg.fun functions are used with gamma, lognormal, reciprocal inverse Gaussian and Gaussian kernel on the motor cycle impact data described in (Silverman 1985). The observations consist of accelerometer reading taken through time in an experimentation on the efficiency of crash helmets. The results in Table 4 agree with the shapes of continuous associated kernels of Part (b) of Figure 1; see also Figure 3. In fact, since the lognormal kernel is well concentrated around the target $x$, it gives the best $R^2$ which is 75.9$\mathrm{\%}$. The gamma and the reciprocal inverse Gaussian kernels give similar $R^2$ in the order $73\mathrm{\%}$. Although the Gaussian kernel is well concentrated on the target, it gives the lower result of $R^2=70.90\mathrm{\%}$. This is mainly due to the symmetry of the kernel which cannot change its shapes according to the target.

Table 4: Some expected values of R${}^2$ of nonparametric regressions of the motor cycle impact data (Silverman 1985) by some continuous associated kernels.

	Gamma	Lognormal	Rec. Inv. Gaussian	Gaussian
R${}^2$	$0.7320$	$0.7591$	$0.7328$	$0.7090$

Figure 3: Nonparametric regressions of the motors cycle impact data (Silverman 1985) by some continuous associated kernels.

Bayesian global for binomial kernel

Using Bayes theorem, the joint posterior distribution of $h$ given the observations is $$\pi (h|X_1,X_2,\dots ,X_n)\propto h^{\alpha -1}(1-h{)}^{\beta -1}{\left(\frac{1}{2}{\displaystyle \sum _{i=1}^n\{y_i-{\hat{m}}_{-i}(X_i){\}}^2+b}\right)}^{-(n+2a)/2},$$ where $\propto$ denotes proportional, the reals $a$ and $b$ are the parameters of the inverse gamma distribution IG$(a,b)$, and $\alpha$ and $\beta$ those of the beta distribution $\mathcal{B}e(\alpha ,\beta )$. The estimate ${\hat{h}}_{bay}$ of the smoothing parameter $h$ is given by Markov chain Monte Carlo (MCMC) techniques with Gibbs sampling: $${\hat{h}}_{bay}=\frac{1}{N-N_0}{\displaystyle \sum _{N_0+1}^N{h}^{(t)},}$$ where $N_0$ is the burn-in period and $N$ the number of iterations; see (Zougab et al. 2014a) for further details. It will be implemented in a future version of the Ake package.

5 Summary and final remarks

The Ake package offers easy tools for R users whose research involves kernel estimation of density functions and/or regression functions through associated kernels that are capable of handling all categorical, count and real positive datasets. Figure 1 shows the importance of the associated kernel choice as well as the bandwidth selection. In fact, symmetric (e.g., Gausssian) kernel estimators (respectively empirical estimators) are not suitable for bounded or positive continuous datasets (respectively discrete small samples). We then need an appropriate associated kernel. The binomial kernel is suitable for small size count data while the discrete triangular or the naive kernel are more indicated for large sample sizes. In continuous cases, the lognormal and gamma kernels give the best estimation for positive data while the extended beta is suitable for any compact support.

This package includes various continuous and discrete associated kernels. It also contains functions to handle the bandwidth selection problems through cross-validation, local and global Bayesian procedures for binomial kernel and also the adaptive Bayesian procedure for the gamma kernel. In general, Bayesian choices of smoothing parameters will be better than their cross-validation counterparts. Future versions of the package will contain Bayesian methods with other associated kernels. Also, these associated kernels are useful for heavy tailed data p.d.f. estimation and can be added later in the package; see, e.g., (Ziane et al. 2015). The case of multivariate data needs to be taken in consideration; see (Kokonendji and Somé 2015) for p.d.f. estimation and (Somé and Kokonendji 2016) for regression. We think that the Ake package can be of interest to nonparametric practitioners of different applied settings.

6 Acknowledgments

We sincerely thank the CRAN editors and two anonymous reviewers for their valuable comments. Part of this work was done while the first author was at “Laboratoire de Mathématiques de Besançon” as a Visiting Scientist, with the support of “Agence Universitaire de la Francophonie” (AUF).