Let us assume that the support of unknown cumulative distribution function (cdf) F is the positive half-line ℝ₊ = (0,∞). To avoid an edge effect when estimating the density function of F it is common to use kernels with the same support as that of the target distribution. Recently, the constructions with asymmetric kernels have been studied for estimating a probability density function (pdf) defined on ℝ₊. Namely, in Chen (2000) and Scaillet (2004) the sequences of gamma kernels, and inverse and reciprocal inverse gaussian kernels have been used, respectively. See also Mnatsakanov and Ruymgaart (2012), where another varying kernel approach is suggested. Their method is based on the sequence of gamma pdfs with varying shapes.

We propose to use a sequence of inverse gamma kernels that represent the δ-sequences in L₂- and L₁-norms, see Lemmas 4.1 and 4.2, respectively. The constructions fα∗ and f̂_α considered in (2.3) and (2.4) (called the varying kernel density estimators (vKDEs)) are different from the traditional kernel density estimator (KDE) (see, for example, Parzen (1962), Silverman (1986), and Scott (1992)). They are also different from the ones proposed in Chen (2000) and Scaillet (2004). In the kernel density estimation the convolution is considered with respect to addition as the group operation on the entire real line ℝ and with a fixed kernel. Our constructions in (2.3) and (2.4) turns out to be of kernel type provided that convolution is considered on the space of a positive half-line (ℝ₊, dH) equipped with multiplication as a group operation, and with the Haar measure dH(t) = dt/t (see, for example, (2.7) below). It is worth mentioning that the estimators proposed by Chen (2000) and by Scaillet (2004) cannot be viewed as convolutions as well as the densities on ℝ₊.

In this paper we investigated the Mean Squared Error (MSE) and Mean Integrated Squared Error (MISE) rates of convergence for proposed estimators fα∗ and f̂_α. Note that the shape of an inverse gamma density varies according to the position of a point x at which the pdf f (x) is estimated. This allows automatic changing the “smoothing” degree around the point x. Another feature of the constructions (2.3) and (2.4) are that they have no boundary effects (see Figs. 1 and 2) and they achieve the optimal rate of convergence for MSE and for MISE within the class of non-negative kernel density estimators. Similar results have been derived in papers: Chen (2000) and Scaillet (2004). There are differences regarding the constants appearing in the first order terms only. It is worth mentioning that in contrast with KDE, the asymptotic variances of fα∗(x) and f̂_α(x) have the same form n-4/5f(x)/(2xπ), as α = n^2/5 (see (3.6) in Section 3), that becomes smaller as x increases. Finally, note that in the case of asymmetric gamma kernels (see Chen (2000)), the corresponding variance has the form n-4/5f(x)/(2xπ). In Mnatsakanov and Ruymgaart (2012), the construction similar to (2.1) has been used, and, as a result, another, the so-called moment-density estimate has been proposed, and its asymptotic properties were studied as well.

The paper is organized as follows. In Section 2 the assumptions and the construction of the vKDE are introduced. In Section 3 the MSE of fα∗ and f̂_α are derived, while in Section 4 the MISE and L₁-consistency of f̂_α are investigated. In Section 5 we conducted the simulation study and compared the performances of the estimators f̂_α, fα∗ and the traditional KDE f̂_h.

In this section we outline the main idea that yields vKDEs fα∗ and f̂_α in (2.3) and (2.4), respectively. Assume we would like to recover (approximate) the moment-identifiable distribution F given only the sequence of its moments. About the conditions necessary and sufficient for F to be the moment-identifiable distribution, see, for example, Stoyanov (2000) and references therein. Suppose that all negative order moments of F are finite. Define the operator ℳ by

(ℳF)(j)=∫0∞t-jdF(t)=μj,j=0,1,… and introduce the sequence of operators ℳα-1: (2.1)(Mα-1μ)(x)=1-∑k=0α(αx)kk!∑j=k∞(-αx)j-k(j-k)!μj,x∈ℝ+.

Here μ = {μ_j, j = 0, 1, …} and α → ∞ at a rate to be specified later.

In analysis, the transform (ℳF)(1 − z), where z is a complex variable, is known as the Mellin transform. There is extensive literature investigating the problem of recovering a function from its Mellin transform. See, for instance, Tagliani (2001), Klauder et al. (2001) and Sneddon (1974), among others. In Gzyl and Tagliani (2010), and Mnatsakanov (2008a,b) the problem of recovering the cdf and corresponding density function given the moment sequence of positive orders of underlying distribution has been studied. The investigation of the properties of approximation Mα-1 in (2.1) is beyond the scope of this article and will be conducted in a separate investigation.

To construct the density estimate, at first, let us approximate F by means of Mα-1. A minor modification of an argument in Mnatsakanov and Ruymgaart (2003) yields

Here by →_w we denote the weak convergence of corresponding cdfs.

Now, suppose we are given a sequence X₁, …, X_n of independent and identically distributed positive random variables from the absolutely continuous distribution function F (with pdf f = F′). To estimate F, let us first estimate its negative j-th order moments μ_j, j ≥ 1. Namely, based on (2.2), let us construct the estimate Fα∗ of F by replacing the moment μ_j in (2.1) by its empirical counterpart

Here F̂_n is the empirical cdf of the sample X₁, …, X_n. After a simple algebra, we derive

To compare Fα∗ with the empirical cdf F̂_n, note that Fα∗(x)~F^n(x) as long as α is large. This follows from the fact that for a given X_i and large α: ∑k=0α1k!(αXix)kexp(-αXix)~I{Xi>x}.

Note also that Fα∗(x) is a continuous function of x, hence, to estimate the density f (x) one can take the derivative of Fα∗(x): (2.3)fα∗(x)=1n∑i=1n1Xi·1Γ(α)(αXix)αexp(-αXix), and choose α = α(n) → ∞ as n → ∞. The problem of optimal choice of parameter α will be specified later. Of course fα∗(x)≥0 for each x > 0, and since it is easily seen that ∫0∞fα∗(x)dx=1, the estimator is itself a probability density. The statements similar to the ones obtained in Sections 3 and 4 are valid for fα∗ as well (see for example, Theorem 3.2). To simplify the calculations below and to reduce the bias of fα∗, let us use the modified version of fα∗. Namely, let us increase the shape parameter of the inverse gamma kernel presented in the right hand side of (2.3) by one. Denote Si,x:=1XiLα(xXi), where L_α(u) = (αu)^α+1/Γ (α + 1) exp(−αu), u ∈ ℝ₊, and consider

Throughout the proposed estimator will be considered at a fixed point x > 0, where f (x) > 0. Also, we will assume that F (0) = 0, and the underlying density satisfies

Besides, let us denote by g(·, a_k, b_k) the inverse gamma density with the shape a_k = k(α + 2) − 1 and the rate b_k = k αx parameters, respectively. Namely

The mean ξ_k and variance σk2 of g(·; a_k, b_k) have the following expressions, respectively: ξk=bkak-1=kαxk(α+2)-2,σk2=bk2(ak-1)2·(ak-2)=k2α2x2{k(α+2)-2}2·{k(α+2)-3}.

Note also that the mean of f̂_α(x) can be written as the convolution operator on (ℝ₊, dH): (2.7)fα(x)=Ef^α(x)=∫0∞Lα(x/t)f(t)dH(t),x∈ℝ+, where dH(t) = dt/t. In Lemmas 4.1 and 4.2, see Section 4, it is proved that the sequence of functions {(1/t) L_α(·/t), t ∈ ℝ₊, α ∈ ℕ} with L_α(·) defined in (2.4) forms the δ-sequences in L₁- and L₂-norms, as α → ∞.

Without explicit reference it will be assumed that all the conditions in Section 2 are satisfied. Let us study the bias and the second moment of the estimator f̂_α. We have

In particular, for k = 1: (3.2)Ef^α(x)=ESi,x=∫0∞g(t;a1,b1)f(t)dt=fα(x).

This yields for the bias of f̂_α(x): (3.3)fα(x)-f(x)=Bias{f^α(x)}=∫0∞g(t;a1,b1){f(t)-f(x)}dt=∫0∞g(t;a1,b1){f(x)+(t-x)f′(x)+12(t-x)2f″(t∼)-f(x)}dt=12∫0∞(t-x)2g(t;a1,b1)f″(x)dt+12∫0∞(t-x)2g(t;a1,b1){f″(t∼)-f″(x)}dt=12·x2α-1·f″(x)+o(1α),asn→∞.

For the variance we have

Applying (3.1) for k = 2 and Bα=α-12-(2α+3)Γ(2α+3)[Γ(α+1)]-2~α1/2/(2π), as α → ∞, yields

Inserting (3.3) and (3.5) in (3.4) we obtain

Finally, combining (3.3) and (3.6) leads to the MSE of f̂_α(x): (3.7)MSE{f^α(x)}=Var{f^α(x)}+Bias2{f^α(x)}=α2nπf(x)x+14x4(α-1)2{f″(x)}2+o(αn)+o(1α2).

For optimal rates we may take

By substitution (3.8) into (3.7) we find

Here we have assumed that the pdf f has a continuous and bounded second derivative f″ (condition (2.5)). The following statement is valid.

In this section we study the performances of fα∗ and f̂_α; defined in (2.3) and (2.4), respectively. In particular, we compare them with KDE f̂_h when the kernel function K is assumed to be a standard normal density function. Let us consider the case when the optimal choice of h, h = h_cv, is based on the least-squares cross validation (CV) algorithm that minimizes the expression M₁(h) defined by Eq. (3.39) in Silverman (1986).

In our simulation studies we plotted the curves of vKDEs f̂_α; and fα∗, when the optimal α = α_cv and, respectively, α=αcv∗, are chosen via the least-squares CV algorithm as well (cf. with Mnatsakanov and Ruymgaart (2012)), and compared them with corresponding curve of KDE f̂_h, when h = h_cv (see Figs. 1 and 2). In particular, we simulated the r.v.’s X_i, i = 1, …, n, from two different distributions: Log-normal (0, 1) and Gamma (2, 1) with different sample sizes n = 200k, 1 ≤ k ≤ 4. In addition, we repeated these simulations N = 500 times and studied the performances of f̂_α, fα∗, and f̂_h using the MISE. Namely, we used the estimated MISE: MISE^:=E^(ISE){f^}=1N∑j=1N∫0∞∣f^(j)(x)-f(x)∣2dx.

Here the expectation Ê is calculated with respect to the empirical cdf of N = 500 values of ISEs, while f̂_(j) denotes the vKDEs or KDE used on the j-th replication. The optimal α = α_cv minimizes the expression M₂(α), i.e.

(5.1)αcv=argminαM2(α)=argminα[∫0∞[f^α(x)]2dx-2∫0∞f^α(x)dF^n(x)], where α ∈ {1, …, 40} for each n = 200k, 1 ≤ k ≤ 4. In the second term of the right hand side of (5.1) let us apply the leave-one-out construction instead of f̂_α. This yields the following expression of

In the case of vKDE fα∗, we choose the optimal CV parameter α=αcv∗ that minimizes the function

During the simulation study, we found out that MISE^s of vKDEs are decreasing functions of n when the parameters α = α_cv, α=αcv∗, and α = n^2/5. In Table 1, we recorded the values of α_cv, αcv∗, and h_cv and corresponding MISE^ for Log-normal (0, 1) and Gamma (2, 1) distributions for four different sample sizes. We see that the values of MISE^ for fα∗ are smaller than corresponding values of MISE^ for f̂_α and f̂_h. To illustrate the performances of vKDEs graphically, we plotted the graphs of estimators f̂_αcv (the dashed curves) with α_cv = 11 and 24 and f̂_h (the dotted curve) with h = h_cv, in Fig. 1(a) and 2(a) when the sampled distributions are Log-normal (0, 1) (with n = 200) and Gamma (2, 1) (with n = 800), respectively. For the same samples, in Fig. 1(b) and 2(b) we plotted the graphs of estimators fαcv∗∗ (the dashed curve) and f̂_h when αcv∗=7 and 18 and h = h_cv, respectively. In each model the sampled pdf f (the solid curve) is plotted as well. Based on the records in Table 1, we conclude that the performances of vKDEs are better compared to the one based on KDE f̂_hcv. After conducting many simulations we can say that the asymptotic behavior of fαcv∗∗ and its modified version f̂_αcv are similar to each other, and their performances around the origin and on the right tail are much better than that of KDE f̂_hcv. For the small sample sizes we suggest to use f̂_αcv instead of f̂_hcv and fαcv∗∗.